Browse code

3.99.1: interventions, death with fdf, user variables

ramon diaz-uriarte (at Phelsuma) authored on 25/06/2022 14:24:13
Showing65 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: 3.5.0
5
-Date: 2021-10-06
4
+Version: 3.99.1
5
+Date: 2022-06-25
6 6
 Authors@R: c(
7 7
 	      person("Ramon", "Diaz-Uriarte", role = c("aut", "cre"),	
8 8
  	   		     email = "rdiaz02@gmail.com"),
... ...
@@ -13,6 +13,14 @@ Authors@R: c(
13 13
               person("Mark", "Taylor", role = "ctb",
14 14
 	        email = "ningkiling@gmail.com",
15 15
 		comment = "plot.stream, plot.stacked"),
16
+	      person("Niklas", "Endres", role = "ctb",
17
+	             comment = "vignette examples, freq-dep-fitness time"),
18
+              person("Javier", "Mu~noz Haro", role = "aut",
19
+		comment = "interventions"),
20
+	      person("Alberto", "Gonzalez Klein", role = "aut",
21
+	        comment = "user-specified death rates"),
22
+	      person("Javier", "Lopez Cano", role = "aut",
23
+	        comment = "user-defined variables"),	
16 24
 	      person("Arash", "Partow", role = "ctb", comment = "ExprTk"),
17 25
 	      person("Sophie", "Brouillet", role = "ctb", comment = "MAGELLAN"),
18 26
 	      person("Sebastian", "Matuszewski", role = "ctb",
... ...
@@ -20,6 +28,7 @@ Authors@R: c(
20 28
 	      person("Harry", "Annoni", role = "ctb", comment = "MAGELLAN"),
21 29
 	      person("Luca", "Ferretti", role = "ctb", comment = "MAGELLAN"),	      
22 30
 	      person("Guillaume", "Achaz", role = "ctb", comment = "MAGELLAN"),
31
+	      person("Wolodzko", "Tymoteusz", role = "ctb", comment = "multivariate hypergeometric"),
23 32
 	      person("Guillermo", "Gorines Cordero", role = "ctb",
24 33
 	             comment = "rfitness"),
25 34
 	      person("Ivan", "Lorca Alonso", role = "ctb",
... ...
@@ -90,13 +99,15 @@ Authors@R: c(
90 99
 	      person("Rafael", "Barrero Rodriguez", role = "ctb",
91 100
 	             comment = "vignette examples"),
92 101
 	      person("Silvia", "Talavera Marcos", role = "ctb",
93
-	             comment = "vignette examples"),
94
-	      person("Niklas", "Endres", role = "ctb",
95
-	             comment = "vignette examples, freq-dep-fitness time")
102
+	             comment = "vignette examples")
96 103
 	      )  
97 104
 Author: Ramon Diaz-Uriarte [aut, cre],
98 105
 	Sergio Sanchez-Carrillo [aut],
99 106
 	Juan Antonio Miguel Gonzalez [aut],
107
+	Alberto Gonzalez Klein [aut],
108
+	Javier Mu\~noz Haro [aut],
109
+	Javier Lopez Cano [aut],
110
+	Niklas Endres [ctb],  
100 111
 	Mark Taylor [ctb],
101 112
 	Arash Partow [ctb],
102 113
 	Sophie Brouillet [ctb],
... ...
@@ -104,6 +115,7 @@ Author: Ramon Diaz-Uriarte [aut, cre],
104 115
 	Harry Annoni [ctb],
105 116
 	Luca Ferretti [ctb],
106 117
 	Guillaume Achaz [ctb],
118
+	Tymoteusz Wolodzko [ctb],
107 119
 	Guillermo Gorines Cordero [ctb],
108 120
 	Ivan Lorca Alonso [ctb],
109 121
 	Francisco Mu\~noz Lopez [ctb],
... ...
@@ -140,27 +152,9 @@ Author: Ramon Diaz-Uriarte [aut, cre],
140 152
 	Guillermo Garcia Hoyos [ctb],
141 153
 	Rosalia Palomino Cabrera [ctb],
142 154
 	Rafael Barrero Rodriguez [ctb],
143
-	Silvia Talavera Marcos [ctb],
144
-	Niklas Endres [ctb]  
155
+	Silvia Talavera Marcos [ctb]	
145 156
 Maintainer: Ramon Diaz-Uriarte <rdiaz02@gmail.com>
146
-Description: Functions for forward population genetic simulation in
147
-    asexual populations, with special focus on cancer progression.
148
-    Fitness can be an arbitrary function of genetic interactions between
149
-    multiple genes or modules of genes, including epistasis, order
150
-    restrictions in mutation accumulation, and order effects.  Fitness can
151
-    also be a function of the relative and absolute frequencies of other
152
-    genotypes (i.e., frequency-dependent fitness). Mutation rates can
153
-    differ between genes, and we can include mutator/antimutator genes (to
154
-    model mutator phenotypes). Simulating multi-species scenarios and
155
-    therapeutic interventions is also possible. Simulations use
156
-    continuous-time models and can include driver and passenger genes and
157
-    modules.  Also included are functions for: simulating random DAGs of
158
-    the type found in Oncogenetic Trees, Conjunctive Bayesian Networks,
159
-    and other cancer progression models; plotting and sampling from single
160
-    or multiple realizations of the simulations, including single-cell
161
-    sampling; plotting the parent-child relationships of the clones;
162
-    generating random fitness landscapes (Rough Mount Fuji, House of
163
-    Cards, additive, NK, Ising, and Eggbox models) and plotting them.
157
+Description: Functions for forward population genetic simulation in asexual populations, with special focus on cancer progression. Fitness can be an arbitrary function of genetic interactions between multiple genes or modules of genes, including epistasis, order restrictions in mutation accumulation, and order effects. Fitness (including just birth, just death, or both and death) can also be a function of the relative and absolute frequencies of other genotypes (i.e., frequency-dependent fitness). Mutation rates can differ between genes, and we can include mutator/antimutator genes (to model mutator phenotypes). Simulating multi-species scenarios and therapeutic interventions, including adaptive therapy, is also possible. Simulations use continuous-time models and can include driver and passenger genes and modules. Also included are functions for: simulating random DAGs of the type found in Oncogenetic Trees, Conjunctive Bayesian Networks, and other cancer progression models; plotting and sampling from single or multiple realizations of the simulations, including single-cell sampling; plotting the parent-child relationships of the clones; generating random fitness landscapes (Rough Mount Fuji, House of Cards, additive, NK, Ising, and Eggbox models) and plotting them.
164 158
 biocViews: BiologicalQuestion, SomaticMutation
165 159
 License: GPL (>= 3)
166 160
 URL: https://github.com/rdiaz02/OncoSimul, https://popmodels.cancercontrol.cancer.gov/gsr/packages/oncosimulr/
... ...
@@ -13,7 +13,10 @@ export("oncoSimulPop", "oncoSimulIndiv", "samplePop",
13 13
        "Magellan_stats",
14 14
        "sampledGenotypes"
15 15
      , "POM", "LOD"
16
-     , "diversityPOM", "diversityLOD"
16
+     , "diversityPOM", "diversityLOD",
17
+     "createInterventions",
18
+     "createUserVars",
19
+     "createRules"
17 20
 ##     , "meanCompositionPop", "compositionPop2"
18 21
        )
19 22
 
... ...
@@ -64,7 +64,10 @@ oncoSimulSample <- function(Nindiv,
64 64
                             fixation = NULL,
65 65
                             verbosity  = 1,
66 66
                             showProgress = FALSE,
67
-                            seed = "auto"){
67
+                            seed = "auto",
68
+                            interventions = NULL,
69
+                            userVars = NULL,
70
+                            rules = NULL){
68 71
     ## No longer using mclapply, because of the way we use the limit on
69 72
     ## the number of tries.
70 73
     
... ...
@@ -192,7 +195,10 @@ oncoSimulSample <- function(Nindiv,
192 195
                                mutationPropGrowth = mutationPropGrowth,
193 196
                                detectionProb = detectionProb,
194 197
                                AND_DrvProbExit = AND_DrvProbExit,
195
-                               fixation = fixation)        
198
+                               fixation = fixation,
199
+                               interventions = interventions,
200
+                                userVars = userVars,
201
+                                rules = rules)        
196 202
         if(tmp$other$UnrecoverExcept) {
197 203
             return(f.out.unrecover.except(tmp))
198 204
         }
... ...
@@ -383,7 +389,10 @@ oncoSimulPop <- function(Nindiv,
383 389
                          fixation = NULL,
384 390
                          verbosity  = 0,
385 391
                          mc.cores = detectCores(),
386
-                         seed = "auto") {
392
+                         seed = "auto",
393
+                         interventions = NULL,
394
+                         userVars = NULL,
395
+                         rules = NULL) {
387 396
 
388 397
     if(Nindiv < 1)
389 398
         stop("Nindiv must be >= 1")
... ...
@@ -424,7 +433,10 @@ oncoSimulPop <- function(Nindiv,
424 433
                         mutationPropGrowth = mutationPropGrowth,
425 434
                         detectionProb = detectionProb,
426 435
                         AND_DrvProbExit = AND_DrvProbExit,
427
-                        fixation = fixation),
436
+                        fixation = fixation,
437
+                        interventions = interventions,
438
+                        userVars = userVars,
439
+                        rules = rules),
428 440
                     mc.cores = mc.cores)
429 441
     ## mc.allow.recursive = FALSE ## FIXME: remove?
430 442
                     ## done for covr issue
... ...
@@ -469,8 +481,10 @@ oncoSimulIndiv <- function(fp,
469 481
                            initMutant = NULL,
470 482
                            AND_DrvProbExit = FALSE,
471 483
                            fixation = NULL,
472
-                           seed = NULL
473
-                           ) {
484
+                           seed = NULL,
485
+                           interventions = NULL,
486
+                           userVars = NULL,
487
+                           rules = NULL) {
474 488
     call <- match.call()
475 489
     if(all(c(is_null_na(detectionProb),
476 490
              is_null_na(detectionSize),
... ...
@@ -495,6 +509,9 @@ oncoSimulIndiv <- function(fp,
495 509
     }
496 510
     if(!inherits(fp, "fitnessEffects")) 
497 511
         stop("v.1 functionality has been removed. Please use v.2")
512
+    
513
+    if(!inherits(fp, "fitnessEffects_v3")) 
514
+        fp <- convertFitnessEffects(fp)
498 515
 
499 516
     ## legacies from poor name choices
500 517
     typeFitness <- switch(model,
... ...
@@ -504,6 +521,8 @@ oncoSimulIndiv <- function(fp,
504 521
                           "McFL" = "mcfarlandlog",
505 522
                           "McFarlandLogD" = "mcfarlandlogd",
506 523
                           "McFLD" = "mcfarlandlogd",
524
+                          "Arb" = "arbitrary",
525
+                          "Const" = "constant",
507 526
                           stop("No valid value for model")
508 527
                           )
509 528
     if(max(initSize) < 1)
... ...
@@ -515,7 +534,21 @@ oncoSimulIndiv <- function(fp,
515 534
     }       ##  if ( !(model %in% c("McFL", "McFarlandLog") )) {
516 535
             ## K <- 1 ## K is ONLY used for McFarland; set it to 1, to avoid
517 536
             ##        ## C++ blowing.
518
-
537
+    
538
+    if("deathSpec" %in% names(fp)) {
539
+        if (fp$deathSpec) {
540
+            if (typeFitness != "arbitrary" && typeFitness != "constant") {
541
+                stop("If death is specified in the fitness effects, use Arb or Const model.")
542
+            }
543
+        }
544
+        
545
+        else {
546
+            if (typeFitness == "arbitrary") {
547
+                stop("To use Arb model specify both birth and death in fitness effects.")
548
+            }
549
+        }
550
+    }
551
+    
519 552
     if(typeFitness == "exp") {
520 553
         death <- 1
521 554
         ## mutationPropGrowth <- 1
... ...
@@ -531,7 +564,7 @@ oncoSimulIndiv <- function(fp,
531 564
     }
532 565
 
533 566
     if(minDetectDrvCloneSz == "auto") {
534
-        if(model %in% c("Bozic", "Exp") )
567
+        if(model %in% c("Bozic", "Exp", "Arb", "Const") )
535 568
             minDetectDrvCloneSz <- 0
536 569
         else if (model %in% c("McFL", "McFarlandLog",
537 570
                               "McFLD", "McFarlandLogD")) {
... ...
@@ -610,6 +643,12 @@ oncoSimulIndiv <- function(fp,
610 643
             if(AND_DrvProbExit)
611 644
                 stop("It makes no sense to pass AND_DrvProbExit and a fixation list.")
612 645
         }
646
+
647
+        #if interventions is null, we create an empty list, cos' it will be easier to handle by C++
648
+        if(is_null_na(interventions)){
649
+            interventions = list()
650
+        }
651
+
613 652
         op <- try(nr_oncoSimul.internal(rFE = fp, 
614 653
                                         birth = birth,
615 654
                                         death = death,  
... ...
@@ -643,7 +682,10 @@ oncoSimulIndiv <- function(fp,
643 682
                                         MMUEF = muEF,
644 683
                                         detectionProb = detectionProb,
645 684
                                         AND_DrvProbExit = AND_DrvProbExit,
646
-                                        fixation = fixation),
685
+                                        fixation = fixation,
686
+                                        interventions = interventions,
687
+                                        userVars = userVars,
688
+                                        rules = rules),
647 689
                   silent = !verbosity)
648 690
         objClass <- c("oncosimul", "oncosimul2")
649 691
     ## }
... ...
@@ -1494,7 +1536,11 @@ get.the.time.for.sample <- function(tmp, timeSample, popSizeSample) {
1494 1536
           } else if (length(candidate.time) == 1) {
1495 1537
                 message("Only one sampled time period with mutants.")
1496 1538
                 the.time <- candidate.time
1497
-            } else {
1539
+          } else {
1540
+              ## If we have drivers at t = 2, t= 4, t= 5, ... but not at t=3
1541
+              ## because at t=3 there were no clones with drivers,
1542
+              ## t = 3 is not among the candidate times.
1543
+              
1498 1544
                   the.time <- sample(candidate.time, 1)
1499 1545
               }
1500 1546
     } else {
... ...
@@ -1,8 +1,8 @@
1 1
 # This file was generated by Rcpp::compileAttributes
2 2
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
3 3
 
4
-nr_BNB_Algo5 <- function(rFE, mu_, death, initSize_, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog, MMUEF, full2mutator_, n2, p2, PDBaseline, cPDetect_i, checkSizePEvery, AND_DrvProbExit, fixation_list) {
5
-    .Call('OncoSimulR_nr_BNB_Algo5', PACKAGE = 'OncoSimulR', rFE, mu_, death, initSize_, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog, MMUEF, full2mutator_, n2, p2, PDBaseline, cPDetect_i, checkSizePEvery, AND_DrvProbExit, fixation_list)
4
+nr_BNB_Algo5 <- function(rFE, mu_, death, initSize_, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog, MMUEF, full2mutator_, n2, p2, PDBaseline, cPDetect_i, checkSizePEvery, AND_DrvProbExit, fixation_list, interventions, userVars, rules) {
5
+    .Call('OncoSimulR_nr_BNB_Algo5', PACKAGE = 'OncoSimulR', rFE, mu_, death, initSize_, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog, MMUEF, full2mutator_, n2, p2, PDBaseline, cPDetect_i, checkSizePEvery, AND_DrvProbExit, fixation_list, interventions, userVars, rules)
6 6
 }
7 7
 
8 8
 evalRGenotype <- function(rG, rFE, spPop, verbose, prodNeg, calledBy_, currentTime) {
... ...
@@ -71,18 +71,21 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,
71 71
     
72 72
   
73 73
     x_from <- y_from <- x_to <- y_to <- Change <- muts <-
74
-        label <- fitness <- Type <- NULL
75
-                
74
+        label <- birth <- Type <- NULL
75
+    
76
+    if ("Fitness" %in% colnames(tfm$afe)) {
77
+      colnames(tfm$afe)[which(colnames(tfm$afe) == "Fitness")] <- "Birth"
78
+    }
76 79
                 
77 80
     dd <- data.frame(muts = mutated,
78
-                     fitness = tfm$afe$Fitness,
81
+                     birth = tfm$afe$Birth,
79 82
                      label = tfm$afe$Genotype,
80 83
                      stringsAsFactors = TRUE)
81 84
     cl <- gaj[vv]
82 85
     sg <- data.frame(x_from = mutated[vv[, 1]],
83
-                     y_from = tfm$afe$Fitness[vv[, 1]],
86
+                     y_from = tfm$afe$Birth[vv[, 1]],
84 87
                      x_to = mutated[vv[, 2]],
85
-                     y_to = tfm$afe$Fitness[vv[, 2]],
88
+                     y_to = tfm$afe$Birth[vv[, 2]],
86 89
                      Change = factor(ifelse(cl == 0, "Neutral",
87 90
                                      ifelse(cl > 0, "Gain", "Loss")),
88 91
                                      levels = c("Gain", "Loss", "Neutral")),
... ...
@@ -91,7 +94,7 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,
91 94
     number_ticks <- function(n) {function(limits) pretty(limits, n)}
92 95
     
93 96
     p <- ggplot() +
94
-        xlab("") + ylab("Fitness") + 
97
+        xlab("") + ylab("Birth") + 
95 98
         theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(),
96 99
               panel.grid.minor.x = element_blank()) +
97 100
         geom_segment(data = sg,
... ...
@@ -121,9 +124,9 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,
121 124
     
122 125
     if(show_labels && use_ggrepel) {
123 126
         p <- p + geom_text_repel(data = dd[-c(minF, maxF), ],
124
-                                 aes(x = muts, y = fitness, label = label)) +
127
+                                 aes(x = muts, y = birth, label = label)) +
125 128
             geom_label_repel(data = ddMM,
126
-                             aes(x = muts, y = fitness, label = label, fill = Type),
129
+                             aes(x = muts, y = birth, label = label, fill = Type),
127 130
                              color = "black",
128 131
                              fontface = "bold")
129 132
             ## geom_label_repel(data = dd[maxF, ], aes(x = muts, y = fitness, label = label),
... ...
@@ -141,10 +144,10 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,
141 144
         }
142 145
             
143 146
         p <- p + geom_text(data = dd[-c(minF, maxF), ],
144
-                           aes(x = muts, y = fitness, label = label),
147
+                           aes(x = muts, y = birth, label = label),
145 148
                            vjust = -0.2, hjust = "inward") +
146 149
             geom_label(data = ddMM,
147
-                       aes(x = muts, y = fitness, label = label, fill = Type),
150
+                       aes(x = muts, y = birth, label = label, fill = Type),
148 151
                        vjust = -0.2, hjust = "inward", color = "black",
149 152
                        fontface = "bold")
150 153
             ## geom_label(data = dd[maxF, ], aes(x = muts, y = fitness, label = label),
... ...
@@ -50,10 +50,27 @@ to_Magellan <- function(x, file,
50 50
 to_Fitness_Matrix <- function(x, max_num_genotypes) {
51 51
     ## A general converter. Ready to be used by plotFitnessLandscape and
52 52
     ## Magellan exporter.
53
-
53
+    
54
+    ## Very bad, but there is no other way
55
+    
56
+    
54 57
     ## FIXME: really, some of this is inefficient. Very. Fix it.
55 58
     if( (inherits(x, "genotype_fitness_matrix")) ||
56 59
         ( (is.matrix(x) || is.data.frame(x)) && (ncol(x) > 2) ) ) {
60
+        
61
+        if("Fitness" %in% colnames(x)) {
62
+            afe <- evalAllGenotypes(allFitnessEffects(
63
+                genotFitness = x, frequencyDependentFitness = FALSE
64
+                ##, epistasis = from_genotype_fitness(x)
65
+            ),
66
+            order = FALSE, addwt = TRUE, max = max_num_genotypes)
67
+        } else {
68
+            afe <- evalAllGenotypes(allFitnessEffects(
69
+                genotFitness = x
70
+                ##, epistasis = from_genotype_fitness(x)
71
+            ),
72
+            order = FALSE, addwt = TRUE, max = max_num_genotypes)
73
+        }
57 74
         ## Why this? We go back and forth twice. We need both things. We
58 75
         ## could construct the afe below by appropriately pasting the
59 76
         ## columns names
... ...
@@ -62,11 +79,7 @@ to_Fitness_Matrix <- function(x, max_num_genotypes) {
62 79
 
63 80
         ## This could use fmatrix_to_afe, above!!!
64 81
         ## Major change as of flfast: no longer using from_genotype_fitness
65
-        afe <- evalAllGenotypes(allFitnessEffects(
66
-            genotFitness = x
67
-            ##, epistasis = from_genotype_fitness(x)
68
-        ),
69
-            order = FALSE, addwt = TRUE, max = max_num_genotypes)
82
+        
70 83
 
71 84
         ## Might not be needed with the proper gfm object (so gmf <- x)
72 85
         ## but is needed if arbitrary matrices.
... ...
@@ -85,10 +98,26 @@ to_Fitness_Matrix <- function(x, max_num_genotypes) {
85 98
         x <- x[, c(1, 2)]
86 99
         if(x[1, "Genotype"] != "WT") {
87 100
             ## Yes, because we expect this present below
88
-            x <- rbind(data.frame(Genotype = "WT",
89
-                                  Fitness = 1,
90
-                                  stringsAsFactors = FALSE),
91
-                       x)
101
+            if ("Death" %in% colnames(x)) {
102
+                x <- rbind(data.frame(Genotype = "WT",
103
+                                      Birth = 1, Death = 1,
104
+                                      stringsAsFactors = FALSE),
105
+                           x)
106
+            } else {
107
+                
108
+                if("Fitness" %in% colnames(x)) {
109
+                    x <- rbind(data.frame(Genotype = "WT",
110
+                                          Fitness = 1,
111
+                                          stringsAsFactors = FALSE),
112
+                               x)
113
+                } else {
114
+                    x <- rbind(data.frame(Genotype = "WT",
115
+                                          Birth = 1,
116
+                                          stringsAsFactors = FALSE),
117
+                               x)
118
+                }
119
+            }
120
+            
92 121
         }
93 122
         afe <- x
94 123
         ## in case we pass an evalAllgenotypesfitandmut
... ...
@@ -113,19 +142,29 @@ to_Fitness_Matrix <- function(x, max_num_genotypes) {
113 142
 ## rfitness, do nothing. Well, it actually does something.
114 143
 
115 144
 to_genotFitness_std <- function(x,
116
-                                frequencyDependentFitness = FALSE,
117
-                                frequencyType = frequencyType,
145
+                                frequencyDependentBirth = FALSE,
146
+                                frequencyDependentDeath = FALSE,
147
+                                frequencyDependentFitness = NULL,
148
+                                frequencyType = NA,
149
+                                deathSpec = FALSE,
118 150
                                 simplify = TRUE,
119
-                                min_filter_fitness = 1e-9,
151
+                                min_filter_birth_death = 1e-9,
120 152
                                 sort_gene_names = TRUE) {
121 153
     ## Would break with output from allFitnessEffects and
122 154
     ## output from allGenotypeAndMut
123 155
 
124 156
     if(! (inherits(x, "matrix") || inherits(x, "data.frame")) )
125 157
         stop("Input must inherit from matrix or data.frame.")
126
-
127
-    if(ncol(x) > 2) {
128
-        if (!frequencyDependentFitness){
158
+    
159
+    if (frequencyDependentDeath && !deathSpec) {
160
+        deathSpec = TRUE
161
+        warning("Assumming death is being specified. Setting deathSpec to TRUE.")
162
+    }
163
+    # Has to be 0's and 1's specification without freq_dep_birth and no death
164
+    # or without freq_dep_birth and without freq_dep_death.
165
+    
166
+    if((ncol(x) > 2 && !deathSpec) || ncol(x) > 3) {
167
+        if (!frequencyDependentBirth && !frequencyDependentDeath){
129 168
             if(inherits(x, "matrix")) {
130 169
                 if(!is.numeric(x))
131 170
                     stop("A genotype fitness matrix/data.frame must be numeric.")
... ...
@@ -139,16 +178,57 @@ to_genotFitness_std <- function(x,
139 178
             if(ncol(x) == 0){
140 179
                 stop("You have an empty data.frame")
141 180
             }
142
-            if(!all(unlist(lapply(x[,-ncol(x)], is.numeric)))){
143
-                stop("All columns except last one must be numeric.")
181
+            if (frequencyDependentBirth && frequencyDependentDeath) {
182
+                if(!all(unlist(lapply(x[,-c((ncol(x)-1):ncol(x))], is.numeric)))){
183
+                    stop("All columns except from the last two must be numeric.")
184
+                }
185
+                
186
+                if(is.factor(x[, ncol(x)])) {
187
+                    warning("Last column of genotype birth-death is a factor. ",
188
+                            "Converting to character.")
189
+                    x[, ncol(x)] <- as.character(x[, ncol(x)])
190
+                }
191
+                
192
+                if(is.factor(x[, ncol(x)-1])) {
193
+                    warning("Second last column of genotype birth-death is a factor. ",
194
+                            "Converting to character.")
195
+                    x[, ncol(x)-1] <- as.character(x[, ncol(x)-1])
196
+                }
197
+                if(!all(unlist(lapply(x[, ncol(x)], is.character)))){
198
+                    stop("All elements in last column must be character.")
199
+                }
200
+                
201
+                if(!all(unlist(lapply(x[, ncol(x)-1], is.character)))){
202
+                    stop("All elements in second last column must be character.")
203
+                }
144 204
             }
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)])
205
+            
206
+            else if(frequencyDependentBirth && deathSpec) {
207
+                if(!all(unlist(lapply(x[,-(ncol(x)-1)], is.numeric)))){
208
+                    stop("All columns except from the second last must be numeric.")
209
+                }
210
+                if(is.factor(x[, ncol(x)-1])) {
211
+                    warning("Second last column of genotype birth-death is a factor. ",
212
+                            "Converting to character.")
213
+                    x[, ncol(x)-1] <- as.character(x[, ncol(x)-1])
214
+                }
215
+                if(!all(unlist(lapply(x[, ncol(x)-1], is.character)))){
216
+                    stop("All elements in second last column must be character.")
217
+                }
149 218
             }
150
-            if(!all(unlist(lapply(x[, ncol(x)], is.character)))){
151
-                stop("All elements in last column must be character.")
219
+            
220
+            else if(frequencyDependentDeath || (frequencyDependentBirth && !deathSpec)) {
221
+                if(!all(unlist(lapply(x[,-ncol(x)], is.numeric)))){
222
+                    stop("All columns except from the last one must be numeric.")
223
+                }
224
+                if(is.factor(x[, ncol(x)])) {
225
+                    warning("Last column of genotype birth-death is a factor. ",
226
+                            "Converting to character.")
227
+                    x[, ncol(x)] <- as.character(x[, ncol(x)])
228
+                }
229
+                if(!all(unlist(lapply(x[, ncol(x)], is.character)))){
230
+                    stop("All elements in last column must be character.")
231
+                }
152 232
             }
153 233
         }
154 234
 
... ...
@@ -158,8 +238,14 @@ to_genotFitness_std <- function(x,
158 238
         ## return(genot_fitness_to_epistasis(x))
159 239
         if(any(duplicated(colnames(x))))
160 240
             stop("duplicated column names")
161
-
162
-        cnfl <- which(colnames(x)[-ncol(x)] == "")
241
+        
242
+        if(deathSpec) {
243
+            cnfl <- which(colnames(x)[-c((ncol(x)-1):ncol(x))] == "")
244
+        }
245
+        else {
246
+            cnfl <- which(colnames(x)[-ncol(x)] == "")
247
+        }
248
+        
163 249
         if(length(cnfl)) {
164 250
             freeletter <- setdiff(LETTERS, colnames(x))[1]
165 251
             if(length(freeletter) == 0) stop("Renaiming failed")
... ...
@@ -168,35 +254,73 @@ to_genotFitness_std <- function(x,
168 254
         }
169 255
         if(!is.null(colnames(x)) && sort_gene_names) {
170 256
             ncx <- ncol(x)
171
-            cnx <- colnames(x)[-ncx]
257
+            if(deathSpec) {
258
+                cnx <- colnames(x)[-c((ncx-1):ncx)]
259
+            }
260
+            else {
261
+                cnx <- colnames(x)[-ncx]
262
+            }
263
+            
172 264
             ocnx <- gtools::mixedorder(cnx)
173 265
             if(!(identical(cnx[ocnx], cnx))) {
174 266
                 message("Sorting gene column names alphabetically")
175
-                x <- cbind(x[, ocnx, drop = FALSE], Fitness = x[, (ncx)])
267
+                
268
+                if(!is.null(frequencyDependentFitness)) {
269
+                    x <- cbind(x[, ocnx, drop = FALSE], Fitness = x[, (ncx)])
270
+                } else {
271
+                    x <- cbind(x[, ocnx, drop = FALSE], Birth = x[, (ncx)])
272
+                }
273
+                
176 274
             }
177 275
         }
178 276
 
179 277
         if(is.null(colnames(x))) {
180
-            ncx <- (ncol(x) - 1)
278
+            if(deathSpec) {
279
+                ncx <- (ncol(x) - 2)
280
+            }
281
+            else {
282
+                ncx <- (ncol(x) - 1)
283
+            }
284
+            
181 285
             message("No column names: assigning gene names from LETTERS")
182 286
             if(ncx > length(LETTERS))
183 287
                 stop("More genes than LETTERS; please give gene names",
184 288
                      " as you see fit.")
185
-            colnames(x) <- c(LETTERS[1:ncx], "Fitness")
289
+            if(deathSpec) {
290
+                colnames(x) <- c(LETTERS[1:ncx], "Birth", "Death")
291
+            }
292
+            
293
+            else {
294
+                
295
+                if (!is.null(frequencyDependentFitness)) {
296
+                    colnames(x) <- c(LETTERS[1:ncx], "Fitness")
297
+                } else {
298
+                    colnames(x) <- c(LETTERS[1:ncx], "Birth")
299
+                }
300
+                    
301
+            }
302
+            
303
+        }
304
+        
305
+        if(deathSpec) {
306
+            if(!all(as.matrix(x[, -c((ncol(x)-1):ncol(x))]) %in% c(0, 1) ))
307
+                stop("First ncol - 2 entries not in {0, 1}.")
308
+        }
309
+        else {
310
+            if(!all(as.matrix(x[, -ncol(x)]) %in% c(0, 1) ))
311
+                stop("First ncol - 1 entries not in {0, 1}.")
186 312
         }
187
-        if(!all(as.matrix(x[, -ncol(x)]) %in% c(0, 1) ))
188
-            stop("First ncol - 1 entries not in {0, 1}.")
189 313
 
190 314
     } else {
191 315
 
192 316
         if(!inherits(x, "data.frame"))
193
-            stop("genotFitness: if two-column must be data frame")
317
+            stop("genotFitness: if genotype is specified, it must be data frame")
194 318
         if(ncol(x) == 0){
195 319
             stop("You have an empty data.frame")
196 320
         }
197 321
         ## Make sure no factors
198 322
         if(is.factor(x[, 1])) {
199
-            warning("First column of genotype fitness is a factor. ",
323
+            warning("First column of genotype birth-death is a factor. ",
200 324
                     "Converting to character.")
201 325
             x[, 1] <- as.character(x[, 1])
202 326
         }
... ...
@@ -221,61 +345,147 @@ to_genotFitness_std <- function(x,
221 345
         if( emarker || ( (!omarker) && (!emarker) && (!nogoodepi)) ) {
222 346
             ## the second case above corresponds to passing just single letter genotypes
223 347
             ## 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")
348
+            if(deathSpec) {
349
+                x <- x[, c(1, 2, 3), drop = FALSE]
350
+                if(!all(colnames(x) == c("Genotype", "Birth", "Death"))) {
351
+                    message("Column names of object not Genotype, Birth and Death.",
352
+                            " Renaming them assuming that is what you wanted")
353
+                    colnames(x) <- c("Genotype", "Birth", "Death")
354
+                }
355
+            }
356
+            else {
357
+                x <- x[, c(1, 2), drop = FALSE]
358
+                if (!is.null(frequencyDependentFitness)) {
359
+                    if(!all(colnames(x) == c("Genotype", "Fitness"))) {
360
+                        message("Column names of object not Genotype and Birth",
361
+                                " Renaming them assuming that is what you wanted")
362
+                        colnames(x) <- c("Genotype", "Fitness")
363
+                    }
364
+                    
365
+                } else {
366
+                    if(!all(colnames(x) == c("Genotype", "Birth"))) {
367
+                        message("Column names of object not Genotype and Birth",
368
+                                " Renaming them assuming that is what you wanted")
369
+                        colnames(x) <- c("Genotype", "Birth")
370
+                    }
371
+                }
372
+                
229 373
             }
374
+            
375
+            
230 376
             if((!omarker) && (!emarker) && (!nogoodepi)) {
231 377
                 message("All single-gene genotypes as input to to_genotFitness_std")
232 378
             }
233 379
             ## Yes, we need to do this to  scale the fitness and put the "-"
234
-            if(frequencyDependentFitness){
380
+            if(frequencyDependentBirth || frequencyDependentDeath){
235 381
                 anywt <- which(x[, 1] == "WT")
236 382
                 if (length(anywt) > 1){
237
-                    stop("WT should not appear more than once in fitness specification")
383
+                    stop("WT should not appear more than once in birth-death specification")
238 384
                 }
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)])
385
+                
386
+                if (frequencyDependentBirth) {
387
+                    if(is.factor(x[, 2])) {
388
+                        warning("Second column of genotype birth-death is a factor. ",
389
+                                "Converting to character.")
390
+                        x[, 2] <- as.character(x[, 2])
391
+                    }
392
+                }
393
+                if (frequencyDependentDeath) {
394
+                    if(is.factor(x[, 3])) {
395
+                        warning("Third column of genotype birth-death is a factor. ",
396
+                                "Converting to character.")
397
+                        x[, 3] <- as.character(x[, 3])
398
+                    }
243 399
                 }
400
+                
244 401
             }
245 402
 
246
-            x <- allGenotypes_to_matrix(x, frequencyDependentFitness)
403
+            x <- allGenotypes_to_matrix(x, frequencyDependentBirth,
404
+                                        frequencyDependentDeath, 
405
+                                        frequencyDependentFitness, deathSpec)
247 406
         }
248 407
     }
249
-    ## And, yes, scale all fitnesses by that of the WT
408
+    ## And, yes, scale all births and deaths by that of the WT
250 409
 
251
-    if (!frequencyDependentFitness){
252
-        whichroot <- which(rowSums(x[, -ncol(x), drop = FALSE]) == 0)
410
+    if (!frequencyDependentBirth && !frequencyDependentDeath){
411
+        if (deathSpec) {
412
+            whichroot <- which(rowSums(x[, -c((ncol(x)-1):ncol(x)), drop = FALSE]) == 0)
413
+        }
414
+        else {
415
+            whichroot <- which(rowSums(x[, -ncol(x), drop = FALSE]) == 0)
416
+        }
417
+        
253 418
         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
419
+            if(deathSpec) {
420
+                warning("No wildtype in the fitness landscape!!! Adding it with birth and death 1.")
421
+                x <- rbind(c(rep(0, ncol(x) - 2), 1, 1), x)
422
+            }
423
+            else {
424
+                warning("No wildtype in the fitness landscape!!! Adding it with birth 1.")
425
+                x <- rbind(c(rep(0, ncol(x) - 1), 1), x)
426
+            }
427
+            
428
+        } else {
429
+            if(x[whichroot, ncol(x)] != 1) {
430
+                if(deathSpec) {
431
+                    warning("Death of wildtype != 1.",
432
+                            " Dividing all deaths by death of wildtype.")
433
+                }
434
+                else {
435
+                    warning("Birth of wildtype != 1.",
436
+                            " Dividing all births by birth of wildtype.")
437
+                }
438
+                
439
+                vwt <- x[whichroot, ncol(x)]
440
+                x[, ncol(x)] <- x[, ncol(x)]/vwt
441
+            }
442
+            
443
+            if (deathSpec && x[whichroot, ncol(x)-1] != 1) {
444
+                warning("Birth of wildtype != 1.",
445
+                        " Dividing all births by birth of wildtype.")
446
+                
447
+                vwt <- x[whichroot, ncol(x)-1]
448
+                x[, ncol(x)-1] <- x[, ncol(x)-1]/vwt
449
+            }
261 450
         }
451
+        
452
+        if(!is.null(frequencyDependentFitness))
453
+            colnames(x)[which(colnames(x) == "Birth")] <- "Fitness"
262 454
     }
263 455
 
264 456
     if(any(is.na(x)))
265 457
         stop("NAs in fitness matrix")
266
-    if(!frequencyDependentFitness) {
458
+    
459
+    if(!frequencyDependentBirth && !frequencyDependentDeath) {
267 460
         if(is.data.frame(x)) 
268 461
             x <- as.matrix(x)
269 462
         stopifnot(inherits(x, "matrix"))
270
-
271
-        if(simplify) {
272
-            x <- x[x[, ncol(x)] > min_filter_fitness, , drop = FALSE]  
463
+    }
464
+    
465
+    if(simplify) {
466
+        if ((!frequencyDependentBirth && !deathSpec) || (!frequencyDependentDeath && deathSpec)) {
467
+            x <- x[x[, ncol(x)] > min_filter_birth_death, , drop = FALSE]
273 468
         }
274
-        class(x) <- c("matrix", "genotype_fitness_matrix")
275
-    } else { ## frequency-dependent fitness
469
+        
470
+        if (!frequencyDependentBirth && deathSpec) {
471
+            x <- x[x[, ncol(x)-1] > min_filter_birth_death, , drop = FALSE]
472
+        }
473
+        
474
+    }
475
+    
476
+    if (!frequencyDependentBirth && !frequencyDependentDeath)
477
+        class(x) <- c("matrix", "genotype_birth_death_matrix")
478
+    
479
+    if(frequencyDependentBirth) { ## frequency-dependent fitness
276 480
         
277 481
         if(frequencyType == "auto"){
278
-            ch <- paste(as.character(x[, ncol(x)]), collapse = "")
482
+            if (deathSpec) {
483
+                ch <- paste(as.character(x[, ncol(x)-1]), collapse = "")
484
+            }
485
+            else {
486
+                ch <- paste(as.character(x[, ncol(x)]), collapse = "")
487
+            }
488
+            
279 489
             if( grepl("f_", ch, fixed = TRUE) ){
280 490
                 frequencyType = "rel"
281 491
                 pattern <- stringr::regex("f_(\\d*_*)*")
... ...
@@ -291,18 +501,62 @@ to_genotFitness_std <- function(x,
291 501
         } else {
292 502
             pattern <- stringr::regex("f_(\\d*_*)*")
293 503
         }
294
-
295
-        regularExpressionVector <-
504
+        
505
+        if(deathSpec) {
506
+            regularExpressionVectorBirth <-
507
+                unique(unlist(lapply(x[, ncol(x)-1],
508
+                                     function(z) {stringr::str_extract_all(string = z,
509
+                                                                           pattern = pattern)})))
510
+            
511
+            if((!all(regularExpressionVectorBirth %in% fVariablesN(ncol(x) - 2, frequencyType))) |
512
+               !(length(intersect(regularExpressionVectorBirth,
513
+                                  fVariablesN(ncol(x) - 2, frequencyType)) >= 1) )){
514
+                stop("There are some errors in birth column")
515
+            }
516
+        }
517
+        else {
518
+            regularExpressionVectorBirth <-
519
+                unique(unlist(lapply(x[, ncol(x)],
520
+                                     function(z) {stringr::str_extract_all(string = z,
521
+                                                                           pattern = pattern)})))
522
+            if((!all(regularExpressionVectorBirth %in% fVariablesN(ncol(x) - 1, frequencyType))) |
523
+               !(length(intersect(regularExpressionVectorBirth,
524
+                                  fVariablesN(ncol(x) - 1, frequencyType)) >= 1) )){
525
+                stop("There are some errors in birth column")
526
+            }
527
+        }
528
+    }
529
+    
530
+    if(frequencyDependentDeath) { ## frequency-dependent fitness
531
+        
532
+        if(frequencyType == "auto"){
533
+            ch <- paste(as.character(x[, ncol(x)]), collapse = "")
534
+            
535
+            if( grepl("f_", ch, fixed = TRUE) ){
536
+                frequencyType = "rel"
537
+                pattern <- stringr::regex("f_(\\d*_*)*")
538
+                
539
+            } else if ( grepl("n_", ch, fixed = TRUE) ){
540
+                frequencyType = "abs"
541
+                pattern <- stringr::regex("n_(\\d*_*)*")
542
+                
543
+            } else { stop("No pattern found when frequencyTypeDeath set to 'auto'") }
544
+            
545
+        } else if(frequencyType == "abs"){
546
+            pattern <- stringr::regex("n_(\\d*_*)*")
547
+        } else {
548
+            pattern <- stringr::regex("f_(\\d*_*)*")
549
+        }
550
+        
551
+        regularExpressionVectorDeath <-
296 552
             unique(unlist(lapply(x[, ncol(x)],
297 553
                                  function(z) {stringr::str_extract_all(string = z,
298 554
                                                                        pattern = pattern)})))
299
-        
300
-        if((!all(regularExpressionVector %in% fVariablesN(ncol(x) - 1, frequencyType))) |
301
-           !(length(intersect(regularExpressionVector,
555
+        if((!all(regularExpressionVectorDeath %in% fVariablesN(ncol(x) - 1, frequencyType))) |
556
+           !(length(intersect(regularExpressionVectorDeath,
302 557
                               fVariablesN(ncol(x) - 1, frequencyType)) >= 1) )){
303
-            stop("There are some errors in fitness column")
558
+            stop("There are some errors in death column")
304 559
         }
305
-
306 560
     }
307 561
     return(x)
308 562
 }
... ...
@@ -373,28 +627,44 @@ genot_fitness_to_epistasis <- function(x) {
373 627
 
374 628
 
375 629
 allGenotypes_to_matrix <- function(x,
376
-                                   frequencyDependentFitness = FALSE) {
630
+                                   frequencyDependentBirth = FALSE,
631
+                                   frequencyDependentDeath = FALSE,
632
+                                   frequencyDependentFitness = NULL,
633
+                                   deathSpec = FALSE) {
377 634
     ## Makes no sense to allow passing order: the matrix would have
378 635
     ## repeated rows. A > B and B > A both have exactly A and B
379 636
 
380 637
     ## Take output of evalAllGenotypes or identical data frame and return
381 638
     ## a matrix with 0/1 in a column for each gene and a final column of
382 639
     ## Fitness
383
-
640
+    
641
+    if(!is.null(frequencyDependentFitness))
642
+        frequencyDependentBirth <- frequencyDependentFitness
643
+    
384 644
     if (is.factor(x[, 1])) {
385 645
         warning(
386
-            "First column of genotype fitness is a factor. ",
646
+            "First column of genotype birth-death is a factor. ",
387 647
             "Converting to character."
388 648
         )
389 649
         x[, 1] <- as.character(x[, 1])
390 650
     }
391
-    if (frequencyDependentFitness) {
392
-        if (is.factor(x[, ncol(x)])) {
651
+    if (frequencyDependentBirth) {
652
+        if (is.factor(x[, 2])) {
653
+            warning(
654
+                "Second column of genotype birth-death is a factor. ",
655
+                "Converting to character."
656
+            )
657
+            x[, 2] <- as.character(x[, 2])
658
+        }
659
+    }
660
+    
661
+    if (frequencyDependentDeath) {
662
+        if (is.factor(x[, 3])) {
393 663
             warning(
394
-                "Second column of genotype fitness is a factor. ",
664
+                "Third column of genotype birth-death is a factor. ",
395 665
                 "Converting to character."
396 666
             )
397
-            x[, ncol(x)] <- as.character(x[, ncol(x)])
667
+            x[, 3] <- as.character(x[, 3])
398 668
         }
399 669
     }
400 670
 
... ...
@@ -402,15 +672,30 @@ allGenotypes_to_matrix <- function(x,
402 672
     anywt <- which(x[, 1] == "WT")
403 673
     if (length(anywt) > 1) stop("More than 1 WT")
404 674
     if (length(anywt) == 1) {
405
-        fwt <- x[anywt, 2]
675
+        bwt <- x[anywt, 2]
676
+        if (deathSpec) {
677
+            dwt <- x[anywt, 3]
678
+        }
679
+        
406 680
         x <- x[-anywt, ]
407 681
         ## Trivial case of passing just a WT?
408 682
     } else {
409
-        if (!frequencyDependentFitness) {
410
-            warning("No WT genotype. Setting its fitness to 1.")
411
-            fwt <- 1
683
+        if (!frequencyDependentBirth && !frequencyDependentDeath) {
684
+            bwt <- 1
685
+            if(deathSpec) {
686
+                dwt <- 1
687
+                warning("No WT genotype. Setting its birth and death to 1.")
688
+            }
689
+            else {
690
+                warning("No WT genotype. Setting its birth to 1.")
691
+            }
692
+            
693
+            
412 694
         } else {
413
-            fwt <- NA
695
+            bwt <- NA
696
+            
697
+            if(deathSpec) 
698
+                dwt <- NA
414 699
             ##   message("No WT genotype in FDF: setting it to 0.")
415 700
         }
416 701
     }
... ...
@@ -436,23 +721,60 @@ allGenotypes_to_matrix <- function(x,
436 721
     for (i in 1:nrow(m)) {
437 722
         m[i, the_match[[i]]] <- 1
438 723
     }
439
-    m <- cbind(m, x[, 2])
440
-    colnames(m) <- c(all_genes, "Fitness")
441
-    if(!is.na(fwt))
442
-        m <- rbind(c(rep(0, length(all_genes)), fwt), m)
724
+    if(deathSpec) {
725
+        m <- cbind(m, x[, 2], x[, 3])
726
+        colnames(m) <- c(all_genes, "Birth", "Death")
727
+    }
728
+    else {
729
+        m <- cbind(m, x[, 2])
730
+        
731
+        if (!is.null(frequencyDependentFitness)) {
732
+            colnames(m) <- c(all_genes, "Fitness")
733
+        } else {
734
+            colnames(m) <- c(all_genes, "Birth")
735
+        }
736
+    }
737
+    
738
+    
739
+    if(!is.na(bwt))
740
+        if(deathSpec) {
741
+            m <- rbind(c(rep(0, length(all_genes)), bwt, dwt), m)
742
+        }
743
+        else{
744
+            m <- rbind(c(rep(0, length(all_genes)), bwt), m)
745
+        }
746
+        
443 747
 
444
-    if (frequencyDependentFitness) {
748
+    if (frequencyDependentBirth || frequencyDependentDeath) {
445 749
         m <- as.data.frame(m)
446 750
         m[, 1:length(all_genes)] <- apply(
447 751
             m[, 1:length(all_genes), drop = FALSE],
448 752
             2,
449 753
             as.numeric
450 754
         )
451
-        m[, ncol(m)] <- as.character(m[, ncol(m)])
755
+        if(frequencyDependentBirth) {
756
+            m[, length(all_genes)+1] <- as.character(m[, length(all_genes)+1])
757
+        }
758
+        else {
759
+            m[, length(all_genes)+1] <- as.numeric(m[, length(all_genes)+1])
760
+        }
761
+        
762
+        if(frequencyDependentDeath) {
763
+            m[, length(all_genes)+2] <- as.character(m[, length(all_genes)+2])
764
+        }
765
+        else if(deathSpec){
766
+            m[, length(all_genes)+2] <- as.numeric(m[, length(all_genes)+2])
767
+        }
768
+        
452 769
     }
453 770
     ## Ensure sorted
454 771
     ## m <- data.frame(m)
455
-    rs <- rowSums(m[, -ncol(m), drop = FALSE])
772
+    if(deathSpec) {
773
+        rs <- rowSums(m[, -c((ncol(m)-1):ncol(m)), drop = FALSE])
774
+    }
775
+    else {
776
+        rs <- rowSums(m[, -ncol(m), drop = FALSE])
777
+    }
456 778
     m <- m[order(rs), , drop = FALSE]
457 779
     ## m <- m[do.call(order, as.list(cbind(rs, m[, -ncol(m)]))), ]
458 780
     return(m)
459 781
new file mode 100644
... ...
@@ -0,0 +1,185 @@
1
+
2
+## Copyright 2013-2021 Ramon Diaz-Uriarte
3
+
4
+## This program is free software: you can redistribute it and/or modify
5
+## it under the terms of the GNU General Public License as published by
6
+## the Free Software Foundation, either version 3 of the License, or
7
+## (at your option) any later version.
8
+
9
+## This program is distributed in the hope that it will be useful,
10
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
11
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12
+## GNU General Public License for more details.
13
+
14
+## You should have received a copy of the GNU General Public License
15
+## along with this program.  If not, see <http://www.gnu.org/licenses/>.
16
+
17
+# we load some functions that are needed to compile
18
+# source("./R/new-restrict.R")
19
+
20
+# this function create the interventions, verifies its correct  specification and returns those interventions
21
+# so they can be processed correctly by C++
22
+createInterventions <- function(interventions, genotFitness, frequencyType = "auto"){
23
+    if((frequencyType == "abs") || (frequencyType == "auto")){
24
+        return (adapt_interventions_to_cpp(verify_interventions(interventions), frequencyType, genotFitness))
25
+    } else if(frequencyType == "rel"){
26
+        stop("Frequency type as relative is not implemented... yet.")
27
+    } else {
28
+        stop("You have specified the freqType wrong. Review it please.")
29
+    }
30
+}
31
+
32
+# this intervention transforms the genotype specification from the user-specified to the C++ one.
33
+# transformIntervention will give more detail about how this works. This "transformation" is done for 
34
+# trigger and what_happens attributes from the intervention.
35
+adapt_interventions_to_cpp <- function(interventions, frequencyType, genotFitness) {
36
+
37
+    for(i in 1:length(interventions)){
38
+        interventions[[i]]$Trigger <- transform_intervention(as.character(interventions[[i]]$Trigger), genotFitness)
39
+        interventions[[i]]$WhatHappens <- transform_intervention(as.character(interventions[[i]]$WhatHappens), genotFitness)
40
+    }
41
+    
42
+    return(interventions)
43
+}
44
+
45
+# this function is the one in charge to transform the genotypes. The user will specify genotypes as A, B, "A,B,C"... etc.
46
+# but in the C++ side, that processes operations, there is no "A", because when fitness is specified, there is an that changes n_A
47
+# to n_1 or n_A_B_C to n_1_2_3. The function that does those "transformations", is all_orders_fv, we just reverse engineered this function
48
+# and borrowed some funcionality so when the user specifies interventions that involve genotypes fitness's, the transformation can be made
49
+transform_intervention <- function(sentence, genotInfo){
50
+    
51
+    prefix <- "n_"
52
+    prefixre <- "^n_"
53
+
54
+    # use functions that handle the change of identifiers, from A -> 1, B->2, A,B -> 1_2
55
+    # functionality "borrowed" by all_orders_fv
56
+    x <- genotInfo$fitnessLandscape[, -ncol(genotInfo$fitnessLandscape), drop = FALSE]
57
+    pasted <- apply(x, 1, function(z) paste(sort(which(z == 1)), collapse = "_"))
58
+    npasted <- apply(x, 1, function(z) paste(sort(colnames(x)[which(z == 1)]), collapse = "_"))
59
+    flvars <- paste0(prefix, pasted)
60
+    names(flvars) <- npasted
61
+    flvars2 <- flvars
62
+    
63
+    # now we have mapped all the original genotypes names with the id'd ones 
64
+    names(flvars2) <- paste0(prefix, names(flvars))
65
+    # we obtain all the posible combinations of the genotypes
66
+    rflvars2 <- rev(flvars2)
67
+    full_rflvars <- all_orders_fv(rflvars2, prefix, prefixre)
68
+    
69
+    # to finish, we replace the previous names with the new ones
70
+    return(stringr::str_replace_all(sentence,
71
+                                    stringr::fixed(full_rflvars)))
72
+}
73
+
74
+## this function checks for inconsistencies in the specification of interventions
75
+verify_interventions <- function(interventions){
76
+
77
+    ## we check if there are interventions with the same ID.
78
+    check_double_id(interventions)
79
+    
80
+    for(i in 1:length(interventions)){
81
+        ## check that the interventions are lists
82
+        print(paste0("Checking intervention: ", interventions[[i]]["ID"]))
83
+        if(is.list(interventions[[i]]) == FALSE){
84
+            stop("Type should be a list of lists")
85
+        }
86
+        
87
+        ## check that exists ID, Trigger, WhatHappens and TimeSensitive in the list 
88
+        if(!exists("ID", interventions[[i]])){
89
+            stop("Attribute ID was not specified.")
90
+        } else if(!exists("Trigger", interventions[[i]])){
91
+            stop("Attribute Trigger was not specified.")
92
+        } else if(!exists("WhatHappens", interventions[[i]])){
93
+            stop("Attribute WhatHappens was not specified.")
94
+        } else if(!exists("Repetitions", interventions[[i]]) || !exists("Periodicity", interventions[[i]]) ){
95
+            stop("Either repetitions or Periodicity must be specified for the intervention")
96
+        }
97
+        
98
+        ## check the data-type of the interventions[[i]]s
99
+        if(!is.character(interventions[[i]]$Trigger) || !is.character(interventions[[i]]$WhatHappens) || !is.character(interventions[[i]]$ID)){
100
+            stop("Interventions are wrongfully specified. Exit...")
101
+        }
102
+
103
+         ## case where user specifies negative periodicity of repetitions
104
+        if (!exists("Periodicity", interventions[[i]]) && interventions[[i]]$Periodicity < 0){
105
+            stop("The periodicity is negative in the intervention")
106
+        }
107
+        if(!exists("Repetitions", interventions[[i]]) && interventions[[i]]$Repetitions < 0){
108
+            stop("The repetitions are negative in the intervention")
109
+        }
110
+        
111
+        ## contemplate the case where only repetitions or only periodicity is specified.
112
+        if(exists("Repetitions", interventions[[i]]) && !(exists("Periodicity", interventions[[i]])) ){
113
+            interventions[[i]]$Periodicity = -1
114
+        } else if(exists("Periodicity", interventions[[i]]) && !(exists("Repetitions", interventions[[i]])) ){ # If the user specifies the periodicity instead of the reps, then TimeSensitive should be setted to "Yes"
115
+            interventions[[i]]$Repetitions = -0.5
116
+        }
117
+
118
+        ## in case user specifies periodicity or reps to infinity, we "cast" it to INT_MAX so C++ understands
119
+        if (exists("Periodicity", interventions[[i]]) && interventions[[i]]$Periodicity == Inf){
120
+            interventions[[i]]$Periodicity <- -1
121
+        }
122
+        if(exists("Repetitions", interventions[[i]]) && interventions[[i]]$Repetitions == Inf){
123
+            interventions[[i]]$Repetitions <- .Machine$integer.max
124
+        }
125
+
126
+        check_what_happens(interventions[[i]]$WhatHappens)
127
+    }
128
+
129
+    return(interventions)
130
+}
131
+
132
+# this function checks that there are no interventions with the same ID.
133
+check_double_id <- function(interventions){
134
+
135
+    if(!is.list(interventions)){
136
+        stop("This should be a list")
137
+    }
138
+
139
+    i <- 1
140
+    buffer <- list()
141
+    while(i <= length(interventions)){
142
+        index_l <- i
143
+        buffer[[i]] <- interventions[[i]]
144
+        j <- 1
145
+        while(j <= length(buffer)){
146
+            if(j!= index_l){
147
+                if(buffer[[j]]$ID == interventions[[i]]$ID){
148
+                    stop("Check the interventions, there are 2 or more that have same IDs")
149
+                }
150
+            }
151
+            j <- j + 1
152
+        } 
153
+        i <- i + 1
154
+    }
155
+}
156
+
157
+# check that the what_happens is correctly specified.
158
+check_what_happens <- function(what_happens){
159
+    # what happens has this form:
160
+    # <genotype_to_apply_some_operation> = <some_operation>
161
+    # we need to assure that the left part is just before the "="
162
+
163
+    string1 <- what_happens
164
+    str_split = strsplit(string1, "\\s+")[[1]]
165
+    # now we check that str_split[[2]] is "="
166
+
167
+    if(str_split[[2]] != "="){
168
+        stop("The specification of WhatHappens is wrong.\n It should be: 
169
+        <genotype_to_apply_some_operation or total_population> = <some_operation>\n Exiting.")
170
+    }
171
+
172
+    flag = FALSE
173
+
174
+    for(s in str_split){
175
+        if(s == "="){
176
+            if(flag == TRUE){
177
+                stop("The specification of WhatHappens is wrong.\n It should be: 
178
+        <genotype_to_apply_some_operation or total_population> = <some_operation>\n Exiting.")
179
+            } else if(flag == FALSE){
180
+                flag = TRUE
181
+            }
182
+        }
183
+    }
184
+}
185
+
... ...
@@ -363,48 +363,128 @@ getGeneIDNum <- function(geneModule, geneNoInt, fitnessLandscape_gene_id,
363 363
 }
364 364
 
365 365
 
366
-## genotFitnes and frequency type -> fitnesLanscapeVariables for FDF and
366
+## genotFitness and frequency type -> fitnessLanscapeVariables for FDF and
367 367
 ##          fitness with numbers, not names
368 368
 ##   Done in a single function since both operations make
369 369
 ##   the same assumptions
370
-create_flvars_fitvars <- function(genotFitness, frequencyType) {
371
-    x <- genotFitness[, -ncol(genotFitness), drop = FALSE]
370
+create_flvars_fitvars <- function(genotFitness, frequencyType, frequencyDependentBirth,
371
+                                  frequencyDependentDeath, frequencyDependentFitness,
372
+                                  deathSpec) {
373
+    if(deathSpec) {
374
+      x <- genotFitness[, -c((ncol(genotFitness)-1):ncol(genotFitness)), drop = FALSE]
375
+      
376
+    }
377
+    else {
378
+      x <- genotFitness[, -ncol(genotFitness), drop = FALSE]
379
+    }
380
+  
372 381
     pasted <- apply(x, 1, function(z) paste(sort(which(z == 1)), collapse = "_"))
373 382
     npasted <- apply(x, 1, function(z) paste(sort(colnames(x)[which(z == 1)]), collapse = "_"))
374
-    if(frequencyType == "abs") {
383
+    
384
+    if(frequencyDependentBirth) {
385
+      
386
+      if(frequencyType == "abs") {
387
+          prefix <- "n_"
388
+          prefixre <- "^n_"
389
+      } else {
390
+          prefix <- "f_"
391
+          prefixre <- "^f_"
392
+      }
393
+      flvarsb <- paste0(prefix, pasted)
394
+      names(flvarsb) <- npasted
395
+  
396
+      ## make sure we get f_1_2 and not f_2_1, etc
397
+      flvars2 <- flvarsb
398
+      names(flvars2) <- paste0(prefix, names(flvarsb))
399
+  
400
+      rmwt <- which(flvars2 == prefix)
401
+      if(length(rmwt)) flvars2 <- flvars2[-rmwt] ## rm this.
402
+  
403
+      ## Need to rev the vector, to ensure larger patterns come first
404
+      ## and to place "f_" as last.
405
+      rflvars2 <- rev(flvars2)
406
+      count_seps <- stringr::str_count(rflvars2, stringr::fixed("_"))
407
+      
408
+      if(any(diff(count_seps) > 0)) {
409
+          warning("flvars not ordered?",
410
+                  "Check the conversion of gene names to numbers in birth spec")
411
+          rflvars2 <- rflvars2[order(count_seps, decreasing = TRUE)]
412
+      }
413
+  
414
+      ## Users can pass in many possible orderings. Get all.
415
+      full_rflvars <- all_orders_fv(rflvars2, prefix, prefixre)
416
+      
417
+      if (!is.null(frequencyDependentFitness)) {
418
+        Fitness_as_fvars <- stringr::str_replace_all(genotFitness$Fitness,
419
+                                                   stringr::fixed(full_rflvars))
420
+      }
421
+      else {
422
+        Birth_as_fvars <- stringr::str_replace_all(genotFitness$Birth,
423
+                                                   stringr::fixed(full_rflvars))
424
+      }
425
+      
426
+    }
427
+    
428
+    if(frequencyDependentDeath) {
429
+      
430
+      if(frequencyType == "abs") {
375 431
         prefix <- "n_"
376 432
         prefixre <- "^n_"
377
-    } else {
433
+      } else {
378 434
         prefix <- "f_"
379 435
         prefixre <- "^f_"
380
-    }
381
-    flvars <- paste0(prefix, pasted)
382
-    names(flvars) <- npasted
383
-
384
-    ## make sure we get f_1_2 and not f_2_1, etc
385
-    flvars2 <- flvars
386
-    names(flvars2) <- paste0(prefix, names(flvars))
387
-
388
-    rmwt <- which(flvars2 == prefix)
389
-    if(length(rmwt)) flvars2 <- flvars2[-rmwt] ## rm this.
390
-
391
-    ## Need to rev the vector, to ensure larger patterns come first
392
-    ## and to place "f_" as last.
393
-    rflvars2 <- rev(flvars2)
394
-    count_seps <- stringr::str_count(rflvars2, stringr::fixed("_"))
395
-    
396
-    if(any(diff(count_seps) > 0)) {
436
+      }
437
+      
438
+      flvarsd <- paste0(prefix, pasted)
439
+      names(flvarsd) <- npasted
440
+      
441
+      ## make sure we get f_1_2 and not f_2_1, etc
442
+      flvars2 <- flvarsd
443
+      names(flvars2) <- paste0(prefix, names(flvarsd))
444
+      
445
+      rmwt <- which(flvars2 == prefix)
446
+      if(length(rmwt)) flvars2 <- flvars2[-rmwt] ## rm this.
447
+      
448
+      ## Need to rev the vector, to ensure larger patterns come first
449
+      ## and to place "f_" as last.
450
+      rflvars2 <- rev(flvars2)
451
+      count_seps <- stringr::str_count(rflvars2, stringr::fixed("_"))
452
+      
453
+      if(any(diff(count_seps) > 0)) {
397 454
         warning("flvars not ordered?",
398
-                "Check the conversion of gene names to numbers in fitness spec")
455
+                "Check the conversion of gene names to numbers in death spec")
399 456
         rflvars2 <- rflvars2[order(count_seps, decreasing = TRUE)]
400
-    }
401
-
402
-    ## Users can pass in many possible orderings. Get all.
403
-    full_rflvars <- all_orders_fv(rflvars2, prefix, prefixre)
404
-    Fitness_as_fvars <- stringr::str_replace_all(genotFitness$Fitness,
457
+      }
458
+      
459
+      ## Users can pass in many possible orderings. Get all.
460
+      full_rflvars <- all_orders_fv(rflvars2, prefix, prefixre)
461
+      Death_as_fvars <- stringr::str_replace_all(genotFitness$Death,
405 462
                                                  stringr::fixed(full_rflvars))
406
-    return(list(flvars = flvars,
407
-                Fitness_as_fvars = Fitness_as_fvars))
463
+      
464
+      if(frequencyDependentBirth) {
465
+        return(list(flvarsb = flvarsb,
466
+                    flvarsd = flvarsd,
467
+                    Birth_as_fvars = Birth_as_fvars,
468
+                    Death_as_fvars = Death_as_fvars))
469
+      }
470
+      else {
471
+        return(list(flvarsd = flvarsd,
472
+                    Death_as_fvars = Death_as_fvars))
473
+      }
474
+      
475
+    }
476
+    
477
+    else  {
478
+      if (!is.null(frequencyDependentFitness)) {
479
+        return(list(flvarsb = flvarsb,
480
+                    Fitness_as_fvars = Fitness_as_fvars))
481
+      } else {
482
+        return(list(flvarsb = flvarsb,
483
+                    Birth_as_fvars = Birth_as_fvars))
484
+      }
485
+     
486
+    }
487
+    
408 488
 }
409 489
 
410 490
 
... ...
@@ -470,7 +550,10 @@ allFitnessORMutatorEffects <- function(rT = NULL,
470 550
                                        genotFitness = NULL,
471 551
                                        ## refFE = NULL,
472 552
                                        calledBy = NULL,
473
-                                       frequencyDependentFitness = FALSE,
553
+                                       frequencyDependentBirth = FALSE,
554
+                                       frequencyDependentDeath = FALSE,
555
+                                       frequencyDependentFitness = NULL,
556
+                                       deathSpec = FALSE,
474 557
                                        frequencyType = "freq_dep_not_used"){
475 558
                                        #spPopSizes = NULL) {
476 559
   ## From allFitnessEffects. Generalized so we deal with Fitness
... ...
@@ -506,7 +589,7 @@ allFitnessORMutatorEffects <- function(rT = NULL,
506 589
            "Is this an attempt to subvert the function?")
507 590
   }
508 591
 
509
-  if(!frequencyDependentFitness) {
592
+  if(!frequencyDependentBirth && !frequencyDependentDeath) {
510 593
     rtNames <- NULL
511 594
     epiNames <- NULL
512 595
     orNames <- NULL
... ...
@@ -626,19 +709,50 @@ allFitnessORMutatorEffects <- function(rT = NULL,
626 709
       ## This makes life simpler in C++:
627 710
       ## In the map, the key is the genotype name, as
628 711
       ## cnn <- colnames(genotFitness)[-ncol(genotFitness)]
629
-      cnn <- 1:(ncol(genotFitness) - 1)
630
-      gfn <- apply(genotFitness[, -ncol(genotFitness), drop = FALSE], 1,
631
-                   function(x) paste(cnn[as.logical(x)],
632
-                                     collapse = ", "))
712
+      if(deathSpec) {
713
+        cnn <- 1:(ncol(genotFitness) - 2)
714
+        gfn <- apply(genotFitness[, -c((ncol(genotFitness)-1):ncol(genotFitness)), drop = FALSE], 1,
715
+                     function(x) paste(cnn[as.logical(x)],
716
+                                       collapse = ", "))
717
+      }
718
+      else {
719
+        cnn <- 1:(ncol(genotFitness) - 1)
720
+        gfn <- apply(genotFitness[, -ncol(genotFitness), drop = FALSE], 1,
721
+                     function(x) paste(cnn[as.logical(x)],
722
+                                       collapse = ", "))
723
+      }
724
+      
633 725
       ## rownames(genotFitness) <- gfn
634
-      fitnessLandscape_df <-
635
-        data.frame(Genotype = gfn,
636
-                   Fitness = genotFitness[, ncol(genotFitness)],
637
-                   stringsAsFactors = FALSE)
638
-      fitnessLandscape_gene_id <- data.frame(
639
-        Gene = colnames(genotFitness)[-ncol(genotFitness)],
640
-        GeneNumID = cnn,
641
-        stringsAsFactors = FALSE)
726
+      if(deathSpec) {
727
+        fitnessLandscape_df <-
728
+          data.frame(Genotype = gfn,
729
+                     Birth = genotFitness[, ncol(genotFitness)-1],
730
+                     Death = genotFitness[, ncol(genotFitness)],
731
+                     stringsAsFactors = FALSE)
732
+        fitnessLandscape_gene_id <- data.frame(
733
+          Gene = colnames(genotFitness)[-c((ncol(genotFitness)-1):ncol(genotFitness))],
734
+          GeneNumID = cnn,
735
+          stringsAsFactors = FALSE)
736
+      }
737
+      else {
738
+        
739
+        if(!is.null(frequencyDependentFitness)) {
740
+          fitnessLandscape_df <-
741
+            data.frame(Genotype = gfn,
742
+                       Fitness = genotFitness[, ncol(genotFitness)],
743
+                       stringsAsFactors = FALSE)
744
+        } else {
745
+          fitnessLandscape_df <-
746
+            data.frame(Genotype = gfn,
747
+                       Birth = genotFitness[, ncol(genotFitness)],
748
+                       stringsAsFactors = FALSE)
749
+        }
750
+        
751
+        fitnessLandscape_gene_id <- data.frame(
752
+          Gene = colnames(genotFitness)[-ncol(genotFitness)],
753
+          GeneNumID = cnn,
754
+          stringsAsFactors = FALSE)
755
+      }
642 756
 
643 757
     }
644 758
 
... ...
@@ -678,32 +792,62 @@ allFitnessORMutatorEffects <- function(rT = NULL,
678 792
     if(!keepInput) {
679 793
       rT <- epistasis <- orderEffects <- noIntGenes <- NULL
680 794
     }
681
-
682
-    out <- list(long.rt = long.rt,
683
-                long.epistasis = long.epistasis,
684
-                long.orderEffects = long.orderEffects,
685
-                long.geneNoInt = geneNoInt,
686
-                geneModule = geneModule,
687
-                gMOneToOne = gMOneToOne,
688
-                geneToModule = geneToModule,
689
-                graph = graphE,
690
-                drv = drv,
691
-                rT = rT,
692
-                epistasis = epistasis,
693
-                orderEffects = orderEffects,
694
-                noIntGenes = noIntGenes,
695
-                fitnessLandscape = genotFitness,
696
-                fitnessLandscape_df = fitnessLandscape_df,
697
-                fitnessLandscape_gene_id = fitnessLandscape_gene_id,
698
-                fitnessLandscapeVariables = vector(mode = "character", length = 0L),
699
-                frequencyDependentFitness = frequencyDependentFitness,
700
-                frequencyType = frequencyType)
795
+    
796
+    if(!is.null(frequencyDependentFitness)) {
797
+      out <- list(long.rt = long.rt,
798
+                  long.epistasis = long.epistasis,
799
+                  long.orderEffects = long.orderEffects,
800
+                  long.geneNoInt = geneNoInt,
801
+                  geneModule = geneModule,
802
+                  gMOneToOne = gMOneToOne,
803
+                  geneToModule = geneToModule,
804
+                  graph = graphE,
805
+                  drv = drv,
806
+                  rT = rT,
807
+                  epistasis = epistasis,
808
+                  orderEffects = orderEffects,
809
+                  noIntGenes = noIntGenes,
810
+                  fitnessLandscape = genotFitness,
811
+                  fitnessLandscape_df = fitnessLandscape_df,
812
+                  fitnessLandscape_gene_id = fitnessLandscape_gene_id,
813
+                  fitnessLandscapeVariables = vector(mode = "character", length = 0L),
814
+                  frequencyDependentFitness = frequencyDependentFitness,
815
+                  frequencyType = frequencyType)
816
+    } else {
817
+    
818
+      out <- list(long.rt = long.rt,
819
+                  long.epistasis = long.epistasis,
820
+                  long.orderEffects = long.orderEffects,
821
+                  long.geneNoInt = geneNoInt,
822
+                  geneModule = geneModule,
823
+                  gMOneToOne = gMOneToOne,
824
+                  geneToModule = geneToModule,
825
+                  graph = graphE,
826
+                  drv = drv,
827
+                  rT = rT,
828
+                  epistasis = epistasis,
829
+                  orderEffects = orderEffects,
830
+                  noIntGenes = noIntGenes,
831
+                  fitnessLandscape = genotFitness,
832
+                  fitnessLandscape_df = fitnessLandscape_df,
833
+                  fitnessLandscape_gene_id = fitnessLandscape_gene_id,
834
+                  fitnessLandscapeVariables = vector(mode = "character", length = 0L),
835
+                  frequencyDependentBirth = frequencyDependentBirth,
836
+                  frequencyDependentDeath = frequencyDependentDeath,
837
+                  frequencyType = frequencyType,
838
+                  deathSpec = deathSpec)
839
+    }
701 840
                 #spPopSizes = vector(mode = "integer", length = 0L)
702 841
     
703 842
     if(calledBy == "allFitnessEffects") {
704
-      class(out) <- c("fitnessEffects")
843
+      if (!is.null(frequencyDependentFitness)) {
844
+        class(out) <- c("fitnessEffects")
845
+      } else {
846
+        class(out) <- c("fitnessEffects", "fitnessEffects_v3")
847
+      }
848
+      
705 849
     } else if(calledBy == "allMutatorEffects") {
706
-      class(out) <- c("mutatorEffects")
850
+      class(out) <- c("mutatorEffects", "mutatorEffects_v3")
707 851
     }
708 852
   } else { ## Frequency-dependent fitness
709 853
 
... ...
@@ -713,39 +857,129 @@ allFitnessORMutatorEffects <- function(rT = NULL,
713 857
       #fitnessLandscape_gene_id <- data.frame()
714 858
       stop("You have a null genotFitness in a frequency dependent fitness situation.")
715 859
     } else {
716
-      cnn <- 1:(ncol(genotFitness) - 1)
717
-      gfn <- apply(genotFitness[, -ncol(genotFitness), drop = FALSE], 1,
718
-                   function(x) paste(cnn[as.logical(x)],
719
-                                     collapse = ", "))
720
-      ## rownames(genotFitness) <- gfn
721
-      fitnessLandscape_df <-
722
-        data.frame(Genotype = gfn,
723
-                   Fitness = genotFitness[, ncol(genotFitness)],
724
-                   stringsAsFactors = FALSE)
860
+      if(deathSpec) {
861
+        cnn <- 1:(ncol(genotFitness) - 2)
862
+        gfn <- apply(genotFitness[, -c((ncol(genotFitness)-1):ncol(genotFitness)), drop = FALSE], 1,
863
+                     function(x) paste(cnn[as.logical(x)],
864
+                                       collapse = ", "))
865
+        
866
+        fitnessLandscape_df <-
867
+          data.frame(Genotype = gfn,
868
+                     Birth = genotFitness[, ncol(genotFitness)-1],
869
+                     Death = genotFitness[, ncol(genotFitness)],
870
+                     stringsAsFactors = FALSE)
871
+        fitnessLandscape_gene_id <- data.frame(
872
+          Gene = colnames(genotFitness)[-c((ncol(genotFitness)-1):ncol(genotFitness))],
873
+          GeneNumID = cnn,
874
+          stringsAsFactors = FALSE)
875
+      }
876
+      else {
877
+        cnn <- 1:(ncol(genotFitness) - 1)
878
+        gfn <- apply(genotFitness[, -ncol(genotFitness), drop = FALSE], 1,
879
+                     function(x) paste(cnn[as.logical(x)],
880
+                                       collapse = ", "))
881
+        
882
+        if (!is.null(frequencyDependentFitness)) {
883
+          fitnessLandscape_df <-
884
+            data.frame(Genotype = gfn,
885
+                       Fitness = genotFitness[, ncol(genotFitness)],
886
+                       stringsAsFactors = FALSE)
887
+        } else {
888
+          fitnessLandscape_df <-
889
+            data.frame(Genotype = gfn,
890
+                       Birth = genotFitness[, ncol(genotFitness)],
891
+                       stringsAsFactors = FALSE)
892
+        }
893
+        
894
+        
895
+        fitnessLandscape_gene_id <- data.frame(
896
+          Gene = colnames(genotFitness)[-ncol(genotFitness)],
897
+          GeneNumID = cnn,
898
+          stringsAsFactors = FALSE)
899
+      }
725 900
 
726 901
       attr(fitnessLandscape_df,'row.names') <-
727 902
         as.integer(attr(fitnessLandscape_df,'row.names'))
728 903
 
729
-      fitnessLandscape_gene_id <- data.frame(
730
-        Gene = colnames(genotFitness)[-ncol(genotFitness)],
731
-        GeneNumID = cnn,
732
-        stringsAsFactors = FALSE)
733
-
734 904
       if(frequencyType == "auto"){
735
-        ch <- paste(as.character(fitnessLandscape_df[, ncol(fitnessLandscape_df)]), collapse = "")
736
-        #print(ch)
737
-        if( grepl("f_", ch, fixed = TRUE) ){
738
-          frequencyType = "rel"
739
-        } else{
740
-          frequencyType = "abs"
905
+        
906
+        if(frequencyDependentBirth && frequencyDependentDeath) {
907
+          
908
+          # We have to make sure both are the same frequencyType
909
+          ch <- paste(as.character(fitnessLandscape_df$Birth), collapse = "")
910
+          ch2 <- paste(as.character(fitnessLandscape_df$Death), collapse = "")
911
+          
912
+          if( grepl("f_", ch, fixed = TRUE) && grepl("f_", ch2, fixed = TRUE)){
913
+            frequencyType = "rel"
914
+          } else if (grepl("n_", ch, fixed = TRUE) && grepl("n_", ch2, fixed = TRUE)){
915
+            frequencyType = "abs"
916
+          } else {
917
+            stop("Inconsistent frequencyType between Birth and Death.
918
+                 Both must be of same frequencyType.")
919
+          }
920
+        }
921
+        
922
+        
923
+        else if(frequencyDependentDeath)  {
924
+          ch <- paste(as.character(fitnessLandscape_df$Death), collapse = "")
925
+          if( grepl("f_", ch, fixed = TRUE)){
926
+            frequencyType = "rel"
927
+          } else{
928
+            frequencyType = "abs"
929
+          }
930
+        }
931
+        
932
+        else if(frequencyDependentBirth) {
933
+          if (!is.null(frequencyDependentFitness)) {
934
+            ch <- paste(as.character(fitnessLandscape_df$Fitness), collapse = "")
935
+          } else {
936
+            ch <- paste(as.character(fitnessLandscape_df$Birth), collapse = "")
937
+          }
938
+          
939
+          if( grepl("f_", ch, fixed = TRUE)){
940
+            frequencyType = "rel"
941
+          } else{
942
+            frequencyType = "abs"
943
+          }
741 944
         }
742
-      } else { frequencyType = frequencyType }
945
+        
946
+      } else { frequencyType = frequencyType}
947
+    
743 948
       ## Wrong: assumes all genotypes in fitness landscape
744 949
       ## fitnessLandscapeVariables = fVariablesN(ncol(genotFitness) - 1, frequencyType)
745
-      stopifnot(identical(genotFitness$Fitness, fitnessLandscape_df$Fitness))
746
-      flvars_and_fitvars <- create_flvars_fitvars(genotFitness, frequencyType)
747
-      fitnessLandscapeVariables <- flvars_and_fitvars$flvars
748
-      Fitness_as_fvars <- flvars_and_fitvars$Fitness_as_fvars
950
+      if(!is.null(frequencyDependentFitness)) {
951
+        stopifnot(identical(genotFitness$Fitness, fitnessLandscape_df$Fitness))
952
+      } else {
953
+        stopifnot(identical(genotFitness$Birth, fitnessLandscape_df$Birth))
954
+      }
955
+      
956
+      if(deathSpec) {
957
+        stopifnot(identical(genotFitness$Death, fitnessLandscape_df$Death))
958
+      }
959
+      flvars_and_fitvars <- create_flvars_fitvars(genotFitness, 
960
+                                                  frequencyType,
961
+                                                  frequencyDependentBirth,
962
+                                                  frequencyDependentDeath,
963
+                                                  frequencyDependentFitness,
964
+                                                  deathSpec)
965
+      
966
+      if(frequencyDependentBirth) {
967
+        
968
+        if (!is.null(frequencyDependentFitness)) {
969
+          fitnessLandscapeVariables <- flvars_and_fitvars$flvarsb
970
+          Fitness_as_fvars <- flvars_and_fitvars$Fitness_as_fvars
971
+        } else {
972
+          birthLandscapeVariables <- flvars_and_fitvars$flvarsb
973
+          Birth_as_fvars <- flvars_and_fitvars$Birth_as_fvars
974
+        }
975
+        
976
+      }
977
+      
978
+      
979
+      if(frequencyDependentDeath) {
980
+        deathLandscapeVariables <- flvars_and_fitvars$flvarsd
981
+        Death_as_fvars <- flvars_and_fitvars$Death_as_fvars
982
+      }
749 983
     }
750 984
 
751 985
     if(!is.null(drvNames)) {
... ...
@@ -763,17 +997,86 @@ allFitnessORMutatorEffects <- function(rT = NULL,
763 997
       ## This is what C++ should consume
764 998
      
765 999
       ## This ought to allow to pass fitness spec as letters. Preserve original
766
-      Fitness_original_as_letters <- fitnessLandscape_df$Fitness
767
-      fitnessLandscape_df$Fitness <- Fitness_as_fvars
768
-
1000
+    if(frequencyDependentBirth) {
1001
+      
1002
+      if (!is.null(frequencyDependentFitness)) {
1003
+        Fitness_original_as_letters <- fitnessLandscape_df$Fitness
1004
+        fitnessLandscape_df$Fitness <- Fitness_as_fvars
1005
+      } else {
1006
+        Birth_original_as_letters <- fitnessLandscape_df$Birth
1007
+        fitnessLandscape_df$Birth <- Birth_as_fvars
1008
+      }
1009
+    }
1010
+      
1011
+    if(frequencyDependentDeath) {
1012
+      Death_original_as_letters <- fitnessLandscape_df$Death
1013
+      fitnessLandscape_df$Death <- Death_as_fvars
1014
+    }
1015
+    if(frequencyDependentDeath && frequencyDependentBirth) {
769 1016
       full_FDF_spec <-
770
-          cbind(genotFitness[, -ncol(genotFitness)]
1017
+        cbind(genotFitness[, -c((ncol(genotFitness)-1):ncol(genotFitness))]
771 1018
               , 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
1019
+              , Genotype_as_letters = genotype_letterslabel(genotFitness[,-c((ncol(genotFitness)-1):ncol(genotFitness))])
1020
+              , Genotype_as_fvarsb = birthLandscapeVariables ## used in C++
1021
+              , Genotype_as_fvarsd = deathLandscapeVariables
1022
+              , Birth_as_fvars = Birth_as_fvars
1023
+              , Birth_as_letters = Birth_original_as_letters
1024
+              , Death_as_fvars = Death_as_fvars
1025
+              , Death_as_letters = Death_original_as_letters
1026
+        )
1027
+
1028
+      }
1029
+      
1030
+      else {
1031
+        
1032
+        if(frequencyDependentDeath) {
1033
+          full_FDF_spec <-
1034
+            cbind(genotFitness[, -c((ncol(genotFitness)-1):ncol(genotFitness))]
1035
+                  , Genotype_as_numbers = fitnessLandscape_df$Genotype
1036
+                  , Genotype_as_letters = genotype_letterslabel(genotFitness[, -c((ncol(genotFitness)-1):ncol(genotFitness))])
1037
+                  , Genotype_as_fvarsd = deathLandscapeVariables ## used in C++
1038
+                  , Death_as_fvars = Death_as_fvars
1039
+                  , Death_as_letters = Death_original_as_letters
1040
+            )
1041
+        } else {
1042
+          if(deathSpec) {
1043
+            full_FDF_spec <-
1044
+              cbind(genotFitness[, -c((ncol(genotFitness)-1):ncol(genotFitness))]
1045
+                    , Genotype_as_numbers = fitnessLandscape_df$Genotype
1046
+                    , Genotype_as_letters = genotype_letterslabel(genotFitness[,-c((ncol(genotFitness)-1):ncol(genotFitness))])
1047
+                    , Genotype_as_fvarsb = birthLandscapeVariables ## used in C++
1048
+                    , Birth_as_fvars = Birth_as_fvars
1049
+                    , Birth_as_letters = Birth_original_as_letters
1050
+              )
1051
+          }
1052
+          
1053
+          else {
1054
+            
1055
+            if (!is.null(frequencyDependentFitness)) {
1056
+              full_FDF_spec <-
1057
+                cbind(genotFitness[, -ncol(genotFitness)]
1058
+                      , Genotype_as_numbers = fitnessLandscape_df$Genotype
1059
+                      , Genotype_as_letters = genotype_letterslabel(genotFitness[, -ncol(genotFitness)])
1060
+                      , Genotype_as_fvarsb = fitnessLandscapeVariables ## used in C++
1061
+                      , Fitness_as_fvars = Fitness_as_fvars
1062
+                      , Fitness_as_letters = Fitness_original_as_letters
1063
+                )
1064
+            } else {
1065
+              full_FDF_spec <-
1066
+                cbind(genotFitness[, -ncol(genotFitness)]
1067
+                      , Genotype_as_numbers = fitnessLandscape_df$Genotype
1068
+                      , Genotype_as_letters = genotype_letterslabel(genotFitness[, -ncol(genotFitness)])
1069
+                      , Genotype_as_fvarsb = birthLandscapeVariables ## used in C++
1070
+                      , Birth_as_fvars = Birth_as_fvars
1071
+                      , Birth_as_letters = Birth_original_as_letters
776 1072
                 )
1073
+            }
1074
+            
1075
+          }
1076
+        }
1077
+      }
1078
+      
1079
+      
777 1080
       rownames(full_FDF_spec) <- 1:nrow(full_FDF_spec)
778 1081
       
779 1082
       ## fitnessLanscape and fitnessLandscape_df are now redundant given
... ...
@@ -781,44 +1084,124 @@ allFitnessORMutatorEffects <- function(rT = NULL,
781 1084
       ## single canonical object used.
782 1085
 
783 1086
       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"
1087
+      if(frequencyDependentBirth) {
1088
+        
1089
+        if (!is.null(frequencyDependentFitness)) {
1090
+          rm(fitnessLandscapeVariables)
1091
+          rm(Fitness_as_fvars)
1092
+          rm(Fitness_original_as_letters)
1093
+        } else {
1094
+          rm(birthLandscapeVariables)
1095
+          rm(Birth_as_fvars)
1096
+          rm(Birth_original_as_letters)
1097
+        }
1098
+        
1099
+      }
792 1100
       
793
-      fitnessLandscape_df <- full_FDF_spec[, c("Genotype_as_numbers",
794
-                                               "Fitness_as_fvars")]
795
-      colnames(fitnessLandscape_df) <- c("Genotype", "Fitness")
1101
+      if (frequencyDependentDeath) {
1102
+        
1103
+        rm(deathLandscapeVariables)
1104
+        rm(Death_as_fvars)
1105
+        rm(Death_original_as_letters)
1106
+        
1107
+        if (frequencyDependentBirth) {
1108
+          
1109
+          fitnessLandscape_df <- full_FDF_spec[, c("Genotype_as_numbers",
1110
+                                                   "Birth_as_fvars",
1111
+                                                   "Death_as_fvars")]
1112
+        } else {
1113
+          
1114
+          fitnessLandscape_df <- cbind(full_FDF_spec["Genotype_as_numbers"],
1115
+                                                   genotFitness["Birth"],
1116
+                                       full_FDF_spec["Death_as_fvars"])
1117
+        }
1118
+        
1119
+        colnames(fitnessLandscape_df) <- c("Genotype", "Birth", "Death")
1120
+        
1121
+      } else {
1122
+        
1123
+        if(deathSpec) {
1124
+        
1125
+          fitnessLandscape_df <- cbind(full_FDF_spec[, c("Genotype_as_numbers",
1126
+                                                         "Birth_as_fvars")], 
1127
+                                       genotFitness[, "Death"])
1128
+                                       
1129
+          colnames(fitnessLandscape_df) <- c("Genotype", "Birth", "Death")
1130
+        
1131
+        } else {
1132
+          
1133
+          if (!is.null(frequencyDependentFitness)) {
1134
+            fitnessLandscape_df <- full_FDF_spec[, c("Genotype_as_numbers",
1135
+                                                     "Fitness_as_fvars")]
1136
+            
1137
+            colnames(fitnessLandscape_df) <- c("Genotype", "Fitness")
1138
+          } else {
1139
+            fitnessLandscape_df <- full_FDF_spec[, c("Genotype_as_numbers",
1140
+                                                     "Birth_as_fvars")]
1141
+            
1142
+            colnames(fitnessLandscape_df) <- c("Genotype", "Birth")
1143
+          }
1144
+          
1145
+        }
1146
+      }
796 1147
       
797 1148
       
798
-      out <- list(long.rt = list(),
799
-                long.epistasis = list(),
800
-                long.orderEffects = list(),
801
-                long.geneNoInt = data.frame(),
802
-                geneModule = defaultGeneModuleDF, ##Trick to pass countGenesFe>2,
803
-                gMOneToOne = TRUE,
804
-                geneToModule = c(Root = "Root"),
805
-                graph = list(),
806
-                drv = drv,
807
-                rT = NULL,
808
-                epistasis = NULL,
809
-                orderEffects = NULL,
810
-                noIntGenes = NULL,
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
815
-                frequencyDependentFitness = frequencyDependentFitness,
816
-                frequencyType = frequencyType,
817
-                full_FDF_spec = full_FDF_spec
818
-                #spPopSizes = spPopSizes
819
-              )
1149
+      if(!is.null(frequencyDependentFitness)) {
1150
+        out <- list(long.rt = list(),
1151
+                    long.epistasis = list(),
1152
+                    long.orderEffects = list(),
1153
+                    long.geneNoInt = data.frame(),
1154
+                    geneModule = defaultGeneModuleDF, ##Trick to pass countGenesFe>2,
1155
+                    gMOneToOne = TRUE,
1156
+                    geneToModule = c(Root = "Root"),
1157
+                    graph = list(),
1158
+                    drv = drv,
1159
+                    rT = NULL,
1160
+                    epistasis = NULL,
1161
+                    orderEffects = NULL,
1162
+                    noIntGenes = NULL,
1163
+                    fitnessLandscape = genotFitness, ## redundant
1164
+                    fitnessLandscape_df = fitnessLandscape_df, ## redundant
1165
+                    fitnessLandscape_gene_id = fitnessLandscape_gene_id, 
1166
+                    ## fitnessLandscapeVariables = NULL, ## now part of full_FDF_spec
1167
+                    frequencyDependentFitness = frequencyDependentFitness,
1168
+                    frequencyType = frequencyType,
1169
+                    full_FDF_spec = full_FDF_spec
1170
+                    #spPopSizes = spPopSizes
1171
+        )
1172
+      } else {
1173
+        out <- list(long.rt = list(),
1174
+                    long.epistasis = list(),
1175
+                    long.orderEffects = list(),
1176
+                    long.geneNoInt = data.frame(),
1177
+                    geneModule = defaultGeneModuleDF, ##Trick to pass countGenesFe>2,
1178
+                    gMOneToOne = TRUE,
1179
+                    geneToModule = c(Root = "Root"),
1180
+                    graph = list(),
1181
+                    drv = drv,
1182
+                    rT = NULL,
1183
+                    epistasis = NULL,
1184
+                    orderEffects = NULL,
1185
+                    noIntGenes = NULL,
1186
+                    fitnessLandscape = genotFitness, ## redundant
1187
+                    fitnessLandscape_df = fitnessLandscape_df, ## redundant
1188
+                    fitnessLandscape_gene_id = fitnessLandscape_gene_id, 
1189
+                    ## fitnessLandscapeVariables = NULL, ## now part of full_FDF_spec
1190
+                    frequencyDependentBirth = frequencyDependentBirth,
1191
+                    frequencyDependentDeath = frequencyDependentDeath,
1192
+                    frequencyType = frequencyType,
1193
+                    deathSpec = deathSpec,
1194
+                    full_FDF_spec = full_FDF_spec
1195
+                    #spPopSizes = spPopSizes
1196
+        )
1197
+      }
1198
+    
1199
+    if (!is.null(frequencyDependentFitness)) {
1200
+      class(out) <- c("fitnessEffects")
1201
+    } else {
1202
+      class(out) <- c("fitnessEffects", "fitnessEffects_v3")
1203
+    }
820 1204
 
821
-    class(out) <- c("fitnessEffects")
822 1205
   }
823 1206
   return(out)
824 1207
 }
... ...
@@ -831,65 +1214,98 @@ allFitnessEffects <- function(rT = NULL,
831 1214
                               drvNames = NULL,
832 1215
                               genotFitness = NULL,
833 1216
                               keepInput = TRUE,
834
-                              frequencyDependentFitness = FALSE,
835
-                              frequencyType = NA) {