Browse code

4.1.1: improved and faster intervention tests

ramon diaz-uriarte (at Phelsuma) authored on 24/11/2022 13:24:15
Showing 12 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: 4.1.0
5
-Date: 2022-10-19
4
+Version: 4.1.1
5
+Date: 2022-11-24
6 6
 Authors@R: c(
7 7
 	      person("Ramon", "Diaz-Uriarte", role = c("aut", "cre"),	
8 8
  	   		     email = "rdiaz02@gmail.com"),
... ...
@@ -1,3 +1,6 @@
1
+Changes in version 4.1 (2022-11-24):
2
+	- Faster, cleaner, and better tests of interventions.
3
+
1 4
 Changes in version 4.0 (for BioC 3.16):
2 5
 	- MAJOR change: users can simulate interventions and adaptive therapy,
3 6
 	  as well as define, track, and use in fitness-dependent specifications
... ...
@@ -50,7 +50,8 @@
50 50
 }
51 51
 
52 52
 \details{
53
-    N/A
53
+  See the vignette for details about differences between intervening on
54
+  the total population or over specific genotypes and when do each occur.
54 55
 }
55 56
 
56 57
 \value{
... ...
@@ -137,4 +138,4 @@
137 138
             Periodicity   = Inf
138 139
         )
139 140
     )
140
-}
141 141
\ No newline at end of file
142
+}
... ...
@@ -669,9 +669,13 @@ of the individual done, and the number of attempts and time used.}
669 669
   * ID: The identifier of the intervention, must be unique.
670 670
   * Trigger: The situation in the simulation that triggers/activates 
671 671
              the intervention.
672
-  * What Happens: "What happens" in the simulation. Basically, once the trigger 
673
-                  is satisfied, this defines how the population is affected by 
674
-                  the intervention.
672
+
673
+  * What Happens: "What happens" in the simulation. Basically, once the
674
+               trigger is satisfied, this defines how the population is
675
+               affected by  the intervention. Please see the vignette
676
+               for details about the differences between when
677
+               interventions that affect a single genotype and those that affect the complete        population occur.
678
+  
675 679
   * Periodicity: Defines the periodicity of the intervention.
676 680
   * Repetitions: Defines the maximum repetitions of each intervention in the simulation.
677 681
 }
... ...
@@ -809,6 +813,8 @@ of the individual done, and the number of attempts and time used.}
809 813
   %% : we assess probability of exiting at every sampling time that (the
810 814
   %% closest possible to) a multiple of \code{checkSizePEvery}.
811 815
 
816
+
817
+  
812 818
 }
813 819
 \value{
814 820
   
... ...
@@ -89,7 +89,7 @@ InterventionsInfo createInterventionsInfo(Rcpp::List interventions,
89 89
         iif = addIntervention(iif, iv);
90 90
     }
91 91
 
92
-    // mapping for the genes and their population
92
+    // mapping for the genotypes and their population
93 93
     iif.mapGenoToPop = evalFVars(fitnessEffects, Genotypes, popParams, true);
94 94
 
95 95
     return iif;
... ...
@@ -115,9 +115,7 @@ bool updatePopulations(InterventionsInfo& iif,
115 115
 //void printIntervention(Intervention i);
116 116
 //void printInterventionsInfo(InterventionsInfo iif);
117 117
 // function that applies hypergeometric progressions to the reduction of the population 
118
-void reduceTotalPopulation(InterventionsInfo& iif, 
119
-                           double target, 
120
-                           double totPopSize);
118
+void reduceTotalPopulation(InterventionsInfo& iif, double target);
121 119
 
122 120
 
123 121
 #endif
... ...
@@ -50,7 +50,9 @@ void parseWhatHappens(InterventionsInfo& iif,
50 50
     } 
51 51
     // now we see if the operation affects total population or genotype popolation
52 52
     std::string leftMostWhatHappens = tokens[0];
53
-    leftMostWhatHappens.erase(std::remove(leftMostWhatHappens.begin(), leftMostWhatHappens.end(), ' '), leftMostWhatHappens.end());
53
+    leftMostWhatHappens.erase(std::remove(leftMostWhatHappens.begin(),
54
+					  leftMostWhatHappens.end(), ' '),
55
+			      leftMostWhatHappens.end());
54 56
     if(leftMostWhatHappens.compare("N") == 0)
55 57
         totalPopFlag = true;
56 58
 
... ...
@@ -75,39 +77,40 @@ void parseWhatHappens(InterventionsInfo& iif,
75 77
         std::string errorMessage = "The expression was imposible to parse.";
76 78
         throw std::invalid_argument(errorMessage);
77 79
     } else{
78
-        // value cant have decimals (can't exist a half-cell)
79
-        double res = floor(expression.value());
80
+      // value cant have decimals (you can't have a half-cell)
81
+      double res = floor(expression.value());
82
+      // once the value is calculated, we must assure if the operation is for the total population
83
+      // or for some specific-genotype
84
+      if (totalPopFlag && (res > N)) {
85
+	throw std::runtime_error("You have specified an intervention that is not allowed: setting the target population size larger than the current value.");
86
+      } 
80 87
 
81
-        // once the value is calculated, we must assure if the operation is for the total population
82
-        // or for some specific-genotype
83
-        if (totalPopFlag && (res > N)) {
84
-            // TODO: Throw exception of some kind, this CANNOT happen by any means
85
-            throw std::runtime_error("You have specified an intervention that is not allowed.");
86
-        } 
88
+      // check that user does not create a WhatHappens that creates population: n_A = n_A * 2
89
+      if(!totalPopFlag){
90
+	try { //check with the left part of the equation finding the value for the n_*
91
+	  const double& value = iif.mapGenoToPop.at(leftMostWhatHappens);
92
+	  if(res > value){
93
+	    Rcpp::Rcerr << "In intervention:" << intervention.id <<
94
+	      " with WhatHappens: " << intervention.what_happens <<
95
+	      ". You cannot intervene to generate more population.";
96
+	    throw std::runtime_error("You have specified an intervention that is not allowed: setting the target population size larger than the current value.");
97
+	  }
98
+	}
99
+	catch (const std::out_of_range&) {
100
+	  Rcpp::Rcout << "Key \"" << leftMostWhatHappens.c_str() << "\" not found" << std::endl;
101
+	}
102
+      }
87 103
 
88
-        // check that user does not create a WhatHappens that creates population: n_A = n_A * 2
89
-        if(!totalPopFlag){
90
-            try { //check with the left part of the equation finding the value for the n_*
91
-                const double& value = iif.mapGenoToPop.at(leftMostWhatHappens);
92
-                if(res > value){
93
-                    Rcpp::Rcerr << "In intervention:" << intervention.id << " with WhatHappens: " << intervention.what_happens << ". You cannot intervene to generate more population.";
94
-                    throw std::runtime_error("You have specified an intervention that is not allowed.");
95
-                }
96
-            }
97
-            catch (const std::out_of_range&) {
98
-                Rcpp::Rcout << "Key \"" << leftMostWhatHappens.c_str() << "\" not found" << std::endl;
99
-            }
100
-        }
101
-
102
-        if(totalPopFlag && res == N){ // this case is absurd, but it might happen, we just return.
103
-            return;
104
-        } else if(totalPopFlag && (res < N)){// reduce total amount of population using hipergeometric distribution
105
-            reduceTotalPopulation(iif, res, totPopSize);
106
-        } else { // update new value for genotype
107
-            std::map<std::string, double>::iterator it = iif.mapGenoToPop.find(leftMostWhatHappens); 
108
-            if(it != iif.mapGenoToPop.end()){
109
-                it->second = res; 
110
-            }
111
-        }
104
+      if(totalPopFlag && res == N){ // this case is absurd, but it might happen, we just return.
105
+	return;
106
+      } else if(totalPopFlag && (res < N)){
107
+	// reduce total amount of population using multivar. hipergeometric distribution
108
+	reduceTotalPopulation(iif, res);
109
+      } else { // update new value for genotype
110
+	std::map<std::string, double>::iterator it = iif.mapGenoToPop.find(leftMostWhatHappens); 
111
+	if(it != iif.mapGenoToPop.end()) {
112
+	  it->second = res; 
113
+	}
114
+      }
112 115
     }
113 116
 }
... ...
@@ -5,42 +5,48 @@
5 5
 // nn by default is equal 1
6 6
 // n array con las poblaciones de cada genotipo y su ncols(diferentes tipos de genotipos en la poblaciĆ³n)
7 7
 // target is the target size, to which number the population would get reduced to.
8
-void reduceTotalPopulation(InterventionsInfo& iif, double target, double totPopSize){
8
+void reduceTotalPopulation(InterventionsInfo& iif, double target){
9 9
     
10
-    // first we take all the population from the structure, and we create a vector with its populations
11
-    std::vector<double> populations;
12
-    Rcpp::NumericMatrix rcpp_mhgeo_distribution; 
13
-    Rcpp::NumericMatrix rcpp_populations_matrix;
14
-    Rcpp::NumericVector rcpp_populations;
15
-    Rcpp::NumericVector rcpp_target;
16
-    //auxiliar variable to check the structure (and the atribute mapGenoToPop) is well populated
17
-    double totalPop = 0.0;
18
-    for(auto map : iif.mapGenoToPop){
19
-        totalPop += map.second;
20
-        populations.push_back(map.second);
21
-    }
22
-    // we convert the vector to something R can understand
23
-    rcpp_populations = Rcpp::wrap(populations);
24
-    rcpp_populations.attr("dim") = Dimension(1, populations.size());
10
+  // first we take all the population from the structure, and we create a vector with its populations
11
+  std::vector<double> populations;
12
+  Rcpp::NumericMatrix rcpp_mhgeo_distribution; 
13
+  Rcpp::NumericMatrix rcpp_populations_matrix;
14
+  Rcpp::NumericVector rcpp_populations;
15
+  Rcpp::NumericVector rcpp_target;
16
+  //auxiliar variable to check the structure (and the atribute mapGenoToPop) is well populated
17
+  double totalPop = 0.0;
18
+  for(auto map : iif.mapGenoToPop){
19
+    totalPop += map.second;
20
+    populations.push_back(map.second);
21
+  }
25 22
 
26
-    rcpp_populations_matrix = Rcpp::wrap(rcpp_populations);
27
-    rcpp_target = Rcpp::wrap(target);
23
+  // // this is never triggered as caught before calling this
24
+  // if (totalPop >= target) {
25
+  //   throw std::runtime_error("Intervention: target population size must be strictly less than existing one");
26
+  // }
27
+  
28
+  // we convert the vector to something R can understand
29
+  rcpp_populations = Rcpp::wrap(populations);
30
+  rcpp_populations.attr("dim") = Dimension(1, populations.size());
31
+
32
+  rcpp_populations_matrix = Rcpp::wrap(rcpp_populations);
33
+  rcpp_target = Rcpp::wrap(target);
28 34
     
29
-    // then, we specify the total genotypes of the populations and we obtain a distribution
30
-    rcpp_mhgeo_distribution = my_rmvhyper(1, rcpp_populations_matrix, rcpp_target);
35
+  // then, we specify the total genotypes of the populations and we obtain a distribution
36
+  rcpp_mhgeo_distribution = my_rmvhyper(1, rcpp_populations_matrix, rcpp_target);
31 37
 
32
-    populations = Rcpp::as<std::vector<double>>(rcpp_mhgeo_distribution);
38
+  populations = Rcpp::as<std::vector<double>>(rcpp_mhgeo_distribution);
33 39
 
34
-    //quick check before creating the matrix CANT BE DONE HERE (totPopSize is not updated until crude is filled)
35
-    //if(totPopSize != totalPop){
36
-    //    throw std::runtime_error("TotalPop != totPopSize, exiting...");
37
-    //}
40
+  //quick check before creating the matrix CANT BE DONE HERE (totPopSize is not updated until crude is filled)
41
+  //if(totPopSize != totalPop){
42
+  //    throw std::runtime_error("TotalPop != totPopSize, exiting...");
43
+  //}
38 44
 
39
-    int i=0;
40
-    for(auto &map : iif.mapGenoToPop){
41
-        //check if it goes out of bounds
42
-        map.second = populations[i];
43
-        i++;
44
-    }
45
+  int i=0;
46
+  for(auto &map : iif.mapGenoToPop){
47
+    //check if it goes out of bounds
48
+    map.second = populations[i];
49
+    i++;
50
+  }
45 51
     
46 52
 }
47 53
deleted file mode 100644
... ...
@@ -1,175 +0,0 @@
1
-inittime <- Sys.time()
2
-cat(paste("\n Starting interventions tests", date(), "\n"))
3
-
4
-## FIXME
5
-## These two tests are extremely computationally intensive,
6
-## and I think could be better tested otherwise.
7
-## And this test has many non-idiomatic R constructs
8
-## and it could probably run in a 1/10 of the time
9
-
10
-test_that("1. Drastically reducing a high-fitness genotype population (McFL) | Trigger depends on T and n_*", {
11
-    set.seed(1)
12
-    df3x <- data.frame(Genotype = c("WT", "B", "R"),
13
-                       Fitness = c("1 + (n_ * 0)",
14
-                                    "1.5",
15
-                                    "1.003 + 0.002 * (n_B > 120)"))
16
-
17
-    afd3 <- allFitnessEffects(genotFitness = df3x,
18
-                          frequencyDependentFitness = TRUE,
19
-                          frequencyType = "abs")
20
-
21
-    ## FIXME: why such periodicity?
22
-    interventions <- list(
23
-        list(
24
-            ID          = "intOverBAffectsR",
25
-            Trigger     = "((T >= 20) and (T <= 70)) and (n_B >= 200)",
26
-            WhatHappens = "n_B = 100",
27
-            Repetitions = 50,
28
-            Periodicity = 0.001
29
-        ),
30
-        list(
31
-            ID          = "intOverTotPop",
32
-            Trigger     = "(T >= 80) and (T <= 85) and (n_B >= 40)",
33
-            WhatHappens = "n_B = 20",
34
-            Repetitions = Inf,
35
-            Periodicity = 0.001
36
-        )
37
-    )
38
-
39
-    interventions <- createInterventions(interventions, afd3)
40
-
41
-
42
-    ep2 <- oncoSimulIndiv(
43
-        afd3,
44
-        model = "McFLD",
45
-        mu = 1e-4,
46
-        sampleEvery = 0.001,
47
-        initSize = c(5000, 10, 300),
48
-        initMutant = c("WT", "B", "R"),
49
-        finalTime = 100,
50
-        onlyCancer = FALSE,
51
-        interventions = interventions,
52
-        ## FIXME: this test occasionally fails
53
-        ## as it goes > 200 s in Windows.
54
-        ## This should not be needed
55
-        max.wall.time = 600
56
-    )
57
-
58
-    ## Why the thresholds? 210 here and 40 below.
59
-
60
-    flag <- FALSE
61
-    i <- 20002
62
-    while(i <= 70001){
63
-        if(ep2$pops.by.time[i, 3:3] >= 210){
64
-            flag <- TRUE
65
-        }
66
-        i <- i + 1
67
-    }
68
-
69
-    testthat::expect_equal(flag, FALSE)
70
-
71
-    # then, between the time intervals, T >= 80 and T<=85
72
-    # we control that the B population
73
-    flag <- FALSE
74
-    i <- 80002
75
-
76
-    ## FIXME: why simulate to 100 time units if we only look up
77
-    ## to row 85000?
78
-    while(i <= 85000){
79
-        if(ep2$pops.by.time[i, 3:3] > 40){
80
-            flag <- TRUE
81
-        }
82
-        i <- i + 1
83
-    }
84
-
85
-    testthat::expect_equal(flag, FALSE)
86
-
87
-    # we plot the simulation when no interventions are specified.
88
-    #plot(ep2, show = "genotypes", type = "line")
89
-})
90
-
91
-
92
-
93
-
94
-
95
-test_that("2. Drastically reducing a high-fitness genotype population (Exp) | Trigger depends on n_*", {
96
-    set.seed(1)
97
-    df3x <- data.frame(Genotype = c("WT", "B", "R"),
98
-                       Fitness = c("1 + (n_ * 0)",
99
-                                   "1.5",
100
-                                   "1.003 + 0.002 * (n_B > 120)"))
101
-
102
-    afd3 <- allFitnessEffects(genotFitness = df3x,
103
-                              frequencyDependentFitness = TRUE,
104
-                              frequencyType = "abs")
105
-
106
-    interventions <- list(
107
-        list(
108
-            ID          = "intOverBAffectsR",
109
-            Trigger     = "((T >= 20) and (T <= 70)) and (n_B >= 200)",
110
-            WhatHappens = "n_B = 100",
111
-            Repetitions = 50,
112
-            Periodicity = 0.001
113
-        ),
114
-        list(
115
-            ID          = "intOverTotPop",
116
-            Trigger     = "(T >= 80) and (T <= 85) and (n_B >= 40)",
117
-            WhatHappens = "n_B = 20",
118
-            Repetitions = Inf,
119
-            Periodicity = 0.001
120
-        )
121
-    )
122
-
123
-    interventions <- createInterventions(interventions, afd3)
124
-
125
-    ep2 <- oncoSimulIndiv(
126
-        afd3,
127
-        model = "Exp",
128
-        mu = 1e-4,
129
-        sampleEvery = 0.001,
130
-        initSize = c(5000, 10, 300),
131
-        initMutant = c("WT", "B", "R"),
132
-        finalTime = 100,
133
-        onlyCancer = FALSE,
134
-        interventions = interventions,
135
-        ## FIXME: This huge wall time should not be necessary.
136
-        ## See above; this is slow as hell in Windows.
137
-        max.wall.time = 600)
138
-
139
-    ## In Macs,
140
-    ##   if (ep2$pops.by.time[i, 3:3] >= 210) {
141
-    ##     flag <- TRUE
142
-    ## }`: argument is of length zero
143
-    ## So only run if not on a Mac
144
-    ## FIXME: this is because the above fails hitting wall time
145
-    ##    if (Sys.info()["sysname"] != "Darwin") {
146
-    flag <- FALSE
147
-    i <- 20002
148
-    while(i <= 70001) {
149
-        if(ep2$pops.by.time[i, 3:3] >= 210) {
150
-            flag <- TRUE
151
-        }
152
-        i <- i + 1
153
-    }
154
-    testthat::expect_equal(flag, FALSE)
155
-
156
-
157
-    ## then, between the time intervals, T >= 80 and T<=85
158
-    ## we control that the B population
159
-    flag <- FALSE
160
-    i <- 80002
161
-    while(i <= 85000) {
162
-        if(ep2$pops.by.time[i, 3:3] > 40){
163
-            flag <- TRUE
164
-        }
165
-        i <- i + 1
166
-    }
167
-
168
-    testthat::expect_equal(flag, FALSE)
169
-    ##    }
170
-})
171
-
172
-cat(paste("\n Ending interventions tests", date(), "\n"))
173
-cat(paste("  Took ", round(difftime(Sys.time(), inittime, units = "secs"), 2), "\n\n"))
174
-rm(inittime)
175
-
... ...
@@ -1,8 +1,5 @@
1 1
 inittime <- Sys.time()
2 2
 cat(paste("\n Starting interventions tests", date(), "\n"))
3
-############################################################################################################################
4
-############################################################################################################################
5
-############################################################################################################################
6 3
 
7 4
 test_that("1. A intervention is created correctly",{
8 5
     fa1 <- data.frame(Genotype = c("A", "B"),
... ...
@@ -22,10 +19,11 @@ test_that("1. A intervention is created correctly",{
22 19
         )
23 20
     )
24 21
 
25
-    testthat::expect_output(createInterventions(interventions, afd3), "Checking intervention: intOverA" )
22
+    testthat::expect_output(createInterventions(interventions, afd3),
23
+                            "Checking intervention: intOverA" )
26 24
     interventions <- createInterventions(interventions, afd3)
27 25
 
28
-    # we check that the transformation of the WhatHappens and Trigger atributte is correct
26
+    ## we check that the transformation of the WhatHappens and Trigger atributte is correct
29 27
     testthat::expect_equal(interventions[[1]]$WhatHappens, "n_1 = n_1 * 0.1")
30 28
     testthat::expect_equal(interventions[[1]]$Trigger, "n_2 >= 5")
31 29
 })
... ...
@@ -45,7 +43,8 @@ test_that("2. Two interventions cannot have the same ID (check_double_id)", {
45 43
             Periodicity = Inf
46 44
         )
47 45
     )
48
-    testthat::expect_error(createInterventions(interventions, afd3), "Check the interventions, there are 2 or more that have same IDs")
46
+    testthat::expect_error(createInterventions(interventions, afd3),
47
+                           "Check the interventions, there are 2 or more that have same IDs")
49 48
 })
50 49
 
51 50
 test_that("3. The attribute WhatHappens is correctly specified (check_what_happens)",{
... ...
@@ -60,7 +59,6 @@ test_that("3. The attribute WhatHappens is correctly specified (check_what_happe
60 59
 
61 60
     testthat::expect_error(createInterventions(interventions, afd3),
62 61
                            "The specification of WhatHappens is wrong.\n It should be:")
63
-    ## <genotype_to_apply_some_operation or total_population> = <some_operation>\n Exiting.")
64 62
 
65 63
     interventions <- list(
66 64
         list(ID           = "intOverA",
... ...
@@ -73,7 +71,6 @@ test_that("3. The attribute WhatHappens is correctly specified (check_what_happe
73 71
 
74 72
     testthat::expect_error(createInterventions(interventions, afd3),
75 73
                            "The specification of WhatHappens is wrong.\n It should be:")
76
-    ##    <genotype_to_apply_some_operation or total_population> = <some_operation>\n Exiting.")
77 74
 
78 75
     interventions <- list(
79 76
         list(ID           = "intOverA",
... ...
@@ -86,7 +83,6 @@ test_that("3. The attribute WhatHappens is correctly specified (check_what_happe
86 83
 
87 84
     testthat::expect_error(createInterventions(interventions, afd3),
88 85
                            "The specification of WhatHappens is wrong.\n It should be:")
89
-    ## <genotype_to_apply_some_operation or total_population> = <some_operation>\n Exiting.")
90 86
 })
91 87
 
92 88
 test_that("4. The user cannot create population in an intervention",{
... ...
@@ -131,7 +127,8 @@ test_that("4. The user cannot create population in an intervention",{
131 127
 
132 128
 
133 129
 test_that("5. Drastically reducing A-genotype population (McFL) | Trigger dependending on T", {
134
-
130
+    seed <- round(runif(1, 1, 1e9))
131
+    
135 132
     fa1 <- data.frame(Genotype = c("A", "B"),
136 133
                     Fitness = c("1.001 + (0*n_A)",
137 134
                                 "1.002"))
... ...
@@ -139,48 +136,52 @@ test_that("5. Drastically reducing A-genotype population (McFL) | Trigger depend
139 136
     afd3 <- allFitnessEffects(genotFitness = fa1,
140 137
                           frequencyDependentFitness = TRUE,
141 138
                           frequencyType = "abs")
142
-    # run the simulation without interventions
139
+
140
+    set.seed(seed)
141
+    ## run the simulation without interventions
143 142
     ep1 <- oncoSimulIndiv(
144
-                    afd3,
145
-                    model = "McFL",
146
-                    mu = 1e-4,
147
-                    initSize = c(20000, 20000),
148
-                    initMutant = c("A", "B"),
149
-                    finalTime = 5.2,
150
-                    sampleEvery = 0.001,
151
-                    onlyCancer = FALSE
152
-                    )
143
+        afd3,
144
+        model = "McFL",
145
+        mu = 1e-4,
146
+        initSize = c(20000, 20000),
147
+        initMutant = c("A", "B"),
148
+        finalTime = 5.2,
149
+        sampleEvery = 0.025,
150
+        onlyCancer = FALSE
151
+    )
152
+
153 153
 
154
-    # now we especify intervention to drastically reduce A population
154
+    ## now we especify intervention to drastically reduce A population
155
+
156
+    interv_fract <- runif(1, 0.01, 0.31)
155 157
     interventions <- list(
156 158
     list(ID           = "intOverA",
157
-        Trigger       = "(T >= 5)",
158
-        WhatHappens   = "n_A = n_A * 0.1",
159
-        Repetitions   = 0,
160
-        Periodicity   = Inf
159
+         Trigger       = "(T >= 3)",
160
+         WhatHappens   =  paste0("n_A = n_A * ", interv_fract),
161
+         Repetitions   = 0,
162
+         Periodicity   = Inf
161 163
     ))
162 164
 
163 165
     interventions <- createInterventions(interventions, afd3)
164 166
 
165
-    # run the simulation WITH interventions
167
+    set.seed(seed)
168
+    ## run the simulation WITH interventions
166 169
     ep2 <- oncoSimulIndiv(
167
-                    afd3,
168
-                    model = "McFL",
169
-                    mu = 1e-4,
170
-                    initSize = c(20000, 20000),
171
-                    initMutant = c("A", "B"),
172
-                    sampleEvery = 0.001,
173
-                    finalTime = 5.2,
174
-                    onlyCancer = FALSE,
175
-		            interventions = interventions
176
-                    )
170
+        afd3,
171
+        model = "McFL",
172
+        mu = 1e-4,
173
+        initSize = c(20000, 20000),
174
+        initMutant = c("A", "B"),
175
+        sampleEvery = 0.025,
176
+        finalTime = 5.2,
177
+        onlyCancer = FALSE,
178
+        interventions = interventions
179
+    )
177 180
 
178
-    # test that the value of the population of A is quite lower once the intervention is made
179
-    testthat::expect_gt(ep2$pops.by.time[4995:4995, 2:2], ep2$pops.by.time[5005:5005, 2:2])
181
+    index <- which(ep2$pops.by.time[,1] %in% ep2$other$interventionTimes)
180 182
 
181
-    # since in the first simulation we do not intervene, the population should be greater that
182
-    # when we intervene
183
-    testthat::expect_lt(ep2$pops.by.time[5005:5005, 2:2], ep1$pops.by.time[5005:5005, 2:2])
183
+    expect_equal(ep2$pops.by.time[index, 2],
184
+                 floor(interv_fract * ep1$pops.by.time[index, 2]))
184 185
 
185 186
 })
186 187
 
... ...
@@ -190,469 +191,448 @@ test_that("5. Drastically reducing A-genotype population (McFL) | Trigger depend
190 191
 
191 192
 test_that("6. Drastically reducing A population (Exp) | Trigger dependending on T", {
192 193
 
194
+    seed <- round(runif(1, 1, 1e9))
195
+
193 196
     fa1 <- data.frame(Genotype = c("WT", "A", "B"),
194
-                    Fitness = c("n_*0",
195
-                                "1.5",
196
-                                "1"))
197
+                      Fitness = c("0 * n_", ## we need an expression
198
+                                  "1.7",
199
+                                  "1.1"))
197 200
 
198 201
     afd3 <- allFitnessEffects(genotFitness = fa1,
199
-                          frequencyDependentFitness = TRUE,
200
-                          frequencyType = "abs")
202
+                              frequencyDependentFitness = TRUE,
203
+                              frequencyType = "abs")
201 204
 
205
+    set.seed(seed)
202 206
     ep1 <- oncoSimulIndiv(
203
-                    afd3,
204
-                    model = "Exp",
205
-                    mu = 1e-4,
206
-                    sampleEvery = 0.001,
207
-                    initSize = c(20000, 20000),
208
-                    initMutant = c("A", "B"),
209
-                    finalTime = 5.2,
210
-                    onlyCancer = FALSE
211
-                    )
207
+        afd3,
208
+        model = "Exp",
209
+        mu = 1e-4,
210
+        sampleEvery = 1, # 0.001,
211
+        initSize = c(20000, 20000),
212
+        initMutant = c("A", "B"),
213
+        finalTime = 10, # 5.2,
214
+        onlyCancer = FALSE
215
+    )
216
+
212 217
 
213
-    # now we especify intervention to drastically reduce A population
218
+    ## now we especify intervention to drastically reduce A population
219
+    interv_fract <- runif(1, 0.01, 0.31)
214 220
     interventions <- list(
215 221
         list(ID           = "intOverA",
216
-            Trigger       = "(T >= 5)",
217
-            WhatHappens   = "n_A = n_A * 0.01",
218
-            Repetitions   = 0,
219
-            Periodicity   = Inf
220
-        )
222
+             Trigger       = "(T >= 5)",
223
+             WhatHappens   = paste0("n_A = n_A * ", interv_fract), ## 0.1
224
+             Repetitions   = 0,
225
+             Periodicity   = Inf
226
+             )
221 227
     )
222 228
 
223 229
     interventions <- createInterventions(interventions, afd3)
224 230
 
231
+    set.seed(seed)
225 232
     ep2 <- oncoSimulIndiv(
226
-                    afd3,
227
-                    model = "Exp",
228
-                    mu = 1e-4,
229
-                    sampleEvery = 0.001,
230
-                    initSize = c(20000, 20000),
231
-                    initMutant = c("A", "B"),
232
-                    finalTime = 5.2,
233
-                    onlyCancer = FALSE,
234
-		            interventions = interventions
235
-                    )
233
+        afd3,
234
+        model = "Exp",
235
+        mu = 1e-4,
236
+        sampleEvery = 1, # 0.001,
237
+        initSize = c(20000, 20000),
238
+        initMutant = c("A", "B"),
239
+        finalTime = 10, # 5.2,
240
+        onlyCancer = FALSE,
241
+        interventions = interventions
242
+    )
243
+
244
+    index <- which(ep2$pops.by.time[,1] %in% ep2$other$interventionTimes)
245
+
246
+    
247
+    last <- nrow(ep1$pops.by.time)
248
+    ## when we do not intervene population of A will be bigger than B,
249
+    ## since it has better fitness
250
+    testthat::expect_gt(ep1$pops.by.time[last, 2],
251
+                        ep1$pops.by.time[last, 3])
252
+
253
+    ## once we intervene we test that the value of the population of A
254
+    ## is quite lower once the intervention is made
255
+    testthat::expect_gt(ep2$pops.by.time[index - 1, 2],
256
+                        ep2$pops.by.time[index, 2])
236 257
 
237
-    # when we do not intervene population of A will be bigger than B, since it has better fitness
238
-    testthat::expect_gt(ep1$pops.by.time[5005:5005, 2:2], ep1$pops.by.time[5005:5005, 3:3])
258
+    ## since in the first simulation we do not intervene,
259
+    ## the population should be greater that
260
+    ## when we intervene
261
+    testthat::expect_lt(ep2$pops.by.time[index, 2],
262
+                        ep1$pops.by.time[index, 2])
239 263
 
240
-    # once we intervene we test that the value of the population of A is quite lower once the intervention is made
241
-    testthat::expect_gt(ep2$pops.by.time[4995:4995, 2:2], ep2$pops.by.time[5005:5005, 2:2])
264
+    ## And this is a much stronger test: the population
265
+    ## of the intervened is exactly the intervention fraction
266
+    ## of the non-intervened
267
+
268
+    expect_equal(ep2$pops.by.time[index, 2],
269
+                 floor(interv_fract * ep1$pops.by.time[index, 2]))
242 270
 
243
-    # since in the first simulation we do not intervene, the population should be greater that
244
-    # when we intervene
245
-    testthat::expect_lt(ep2$pops.by.time[5005:5005, 2:2], ep1$pops.by.time[5005:5005, 2:2])
246 271
 })
247 272
 
248 273
 
249 274
 test_that("7. Intervening over total population (McFL) | Trigger depends on T", {
250 275
     gffd3 <- data.frame(Genotype = c("WT", "A", "B"),
251
-                    Fitness = c("1",
252
-                    "1 + 0.25 * (n_B > 0)",
253
-                    ".9 + 0.4 * (n_A > 0)"
254
-                    ))
255
-    afd3 <- allFitnessEffects(genotFitness = gffd3,
256
-                            frequencyDependentFitness = TRUE,
257
-                            frequencyType = "abs")
276
+                        Fitness = c("1",
277
+                                    "1 + 0.2 * (n_B > 0)",
278
+                                    ".9 + 0.4 * (n_A > 0)"
279
+                                    ))
280
+    aafd3 <- allFitnessEffects(genotFitness = gffd3,
281
+                               frequencyDependentFitness = TRUE,
282
+                               frequencyType = "abs")
283
+
284
+
285
+    interv_fract <- round(runif(1, 0.4, 0.8), 2)
258 286
 
259 287
     interventions = list(
260 288
         list(
261 289
             ID            = "intOverTotPop",
262
-            Trigger       = "T > 40",
263
-            WhatHappens   = "N = N * 0.6",
290
+            Trigger       = "T > 5",
291
+            WhatHappens   = paste0("N = N * ", interv_fract),
264 292
             Repetitions   = 2,
265
-            Periodicity   = 20
293
+            Periodicity   = 5
266 294
         )
267 295
     )
268 296
 
269
-    interventions <- createInterventions(interventions, afd3)
270
-
271
-    sfd3 <- oncoSimulIndiv(afd3,
272
-                           model = "McFLD",
297
+    interventions <- createInterventions(interventions, aafd3)
298
+    
299
+    sfd3 <- oncoSimulIndiv(aafd3,
300
+                           model = "McFL",
273 301
                            onlyCancer = FALSE,
274
-                           finalTime = 100,
302
+                           finalTime = 16,
275 303
                            mu = 1e-4,
276
-                           initSize = 5000,
277
-                           sampleEvery = 0.001,
278
-                           interventions = interventions)
279
-
280
-    # we can check genotype by genotype that when an intervention ocurs, their population lowers
281
-
282
-    indexes <- vector()
283
-    # We get the indexes that match the intervention times in pops.by.time
284
-    for(time in sfd3$other$interventionTimes){
285
-        indexes <- append(indexes, which(sfd3$pops.by.time[,1] == time))
286
-    }
287
-
288
-    # For each intervention time (T = 40, 60, 80)
289
-    for(index in indexes){
290
-        line <- sfd3$pops.by.time[index, ]
291
-        prev_line <- sfd3$pops.by.time[index-1, ]
292
-                                        #Total
293
-        total <- line[2] + line[3] + line[4]
294
-        prev_total <- prev_line[2] + prev_line[3] + prev_line[4]
295
-        testthat::expect_gt(total, prev_total*0.6 - 0.2*prev_total)
296
-        testthat::expect_lt(total, prev_total*0.6 + 0.2*prev_total)
297
-            #Genotype WT
298
-        if((prev_line[2] > 0) & (line[2] > 0)){
299
-            testthat::expect_gte(prev_line[2], line[2])
300
-        }
301
-            #Genotype A
302
-        if((prev_line[3] > 0) & (line[3] > 0)){
303
-            testthat::expect_gte(prev_line[3], line[3])
304
-        }
305
-            #Genotype B
306
-        if((prev_line[4] > 0) & (line[4] > 0)){
307
-            testthat::expect_gte(prev_line[4], line[4])
308
-        }
304
+                           initSize = 2e4, 
305
+                           sampleEvery = 0.025, 
306
+                           interventions = interventions,
307
+                           detectionSize = NA
308
+                           )
309
+
310
+    ## it may happen that, in some simulations, the population collapses, in that case,
311
+    ## pops by time is null, and cannot be checked
312
+
313
+    if (!is.null(sfd3$pops.by.time)) {
314
+        indexes <- which(sfd3$pops.by.time[,1] %in% sfd3$other$interventionTimes)
315
+        total_before <- rowSums(sfd3$pops.by.time[indexes - 1, -1])
316
+        total_after <-  rowSums(sfd3$pops.by.time[indexes, -1])
317
+        reduction <- round(total_after/total_before, 2)
318
+        print(reduction)
319
+        expect_equal(reduction, rep(interv_fract, 3))
309 320
     }
321
+    
310 322
 })
311 323
 
312 324
 
313
-
314 325
 test_that("8. Intervening over total population (Exp) | Trigger depends on T", {
315 326
     gffd3 <- data.frame(Genotype = c("WT", "A", "B"),
316
-                    Fitness = c("1",
317
-                    "1 + 0.2 * (n_B > 0)",
318
-                    ".9 + 0.4 * (n_A > 0)"
319
-                    ))
320
-    afd3 <- allFitnessEffects(genotFitness = gffd3,
321
-                            frequencyDependentFitness = TRUE,
322
-                            frequencyType = "abs")
327
+                        Fitness = c("1",
328
+                                    "1 + 0.2 * (n_B > 0)",
329
+                                    ".9 + 0.4 * (n_A > 0)"
330
+                                    ))
331
+    aafd3 <- allFitnessEffects(genotFitness = gffd3,
332
+                               frequencyDependentFitness = TRUE,
333
+                               frequencyType = "abs")
334
+
335
+
336
+    interv_fract <- round(runif(1, 0.4, 0.8), 2)
323 337
 
324 338
     interventions = list(
325 339
         list(
326 340
             ID            = "intOverTotPop",
327 341
             Trigger       = "T > 10",
328
-            WhatHappens   = "N = N * 0.8",
342
+            WhatHappens   = paste0("N = N * ", interv_fract),
329 343
             Repetitions   = 2,
330 344
             Periodicity   = 10
331 345
         )
332 346
     )
333 347
 
334
-    interventions <- createInterventions(interventions, afd3)
335
-
336
-    sfd3 <- oncoSimulIndiv( afd3,
337
-                            model = "Exp",
338
-                            onlyCancer = FALSE,
339
-                            finalTime = 40,
340
-                            mu = 1e-4,
341
-                            initSize = 5000,
342
-                            sampleEvery = 0.001,
343
-                            interventions = interventions)
344
-
345
-    # it may happen that, in some simulations, the population collapses, in that case,
346
-    # pops by time is null, and cannot be checked
347
-
348
-    # we can check genotype by genotype that when an intervention ocurs, their population lowers
349
-    indexes <- vector()
350
-    # We get the indexes that match the intervention times in pops.by.time
351
-    for(time in sfd3$other$interventionTimes){
352
-        indexes <- append(indexes, which(sfd3$pops.by.time[,1] == time))
353
-    }
354
-
355
-    # For each intervention time (T = 10, 20, 30)
356
-    for(index in indexes){
357
-        line <- sfd3$pops.by.time[index, ]
358
-        prev_line <- sfd3$pops.by.time[index-1, ]
359
-                                        #Total
360
-        total <- line[2] + line[3] + line[4]
361
-        prev_total <- prev_line[2] + prev_line[3] + prev_line[4]
362
-        testthat::expect_gt(total, prev_total*0.8 - 0.2*prev_total)
363
-        testthat::expect_lt(total, prev_total*0.8 + 0.2*prev_total)
364
-            #Genotype WT
365
-        if((prev_line[2] > 0) & (line[2] > 0)){
366
-            testthat::expect_gte(prev_line[2], line[2])
367
-        }
368
-            #Genotype A
369
-        if((prev_line[3] > 0) & (line[3] > 0)){
370
-            testthat::expect_gte(prev_line[3], line[3])
371
-        }
372
-            #Genotype B
373
-        if((prev_line[4] > 0) & (line[4] > 0)){
374
-            testthat::expect_gte(prev_line[4], line[4])
375
-        }
348
+    interventions <- createInterventions(interventions, aafd3)
349
+    
350
+    sfd3 <- oncoSimulIndiv(aafd3,
351
+                           model = "Exp",
352
+                           onlyCancer = FALSE,
353
+                           finalTime = 32,
354
+                           mu = 1e-4,
355
+                           initSize = 2e5, ## 20000,
356
+                           sampleEvery = 1, ## 0.001,
357
+                           interventions = interventions,
358
+                           detectionSize = NA
359
+                           )
360
+
361
+    ## it may happen that, in some simulations, the population collapses, in that case,
362
+    ## pops by time is null, and cannot be checked
363
+
364
+    if (!is.null(sfd3$pops.by.time)) {
365
+        indexes <- which(sfd3$pops.by.time[,1] %in% sfd3$other$interventionTimes)
366
+        total_before <- rowSums(sfd3$pops.by.time[indexes - 1, -1])
367
+        total_after <-  rowSums(sfd3$pops.by.time[indexes, -1])
368
+        reduction <- round(total_after/total_before, 2)
369
+        expect_equal(reduction, rep(interv_fract, 3))
376 370
     }
371
+    
377 372
 })
378 373
 
379
-# test 9 and 10 found in test.Z-intervention.R
374
+## test 9 and 10 found in test.Z-intervention.R
380 375
 
381
-test_that("11. Intervening over 4 genotypes both over specific genotype and total population (McFL) | Trigger depends on N",{
376
+test_that("11. Intervening over 4 genotypes both over specific genotype and total population (Exp) | Trigger depends on N",{
382 377
     df3x <- data.frame(Genotype = c("A", "B", "C", "D", "E"),
383
-                      Fitness = c("1",
384
-                                  "1.01 + (0 * n_A)",
385
-                                  "1.1",
386
-                                  "1.09",
387
-                                  "1.07"))
378
+                       Fitness = c("1",
379
+                                   "1.01 + (0 * n_A)",
380
+                                   "1.1",
381
+                                   "1.19",
382
+                                   "1.17"))
388 383
 
389 384
     afd3 <- allFitnessEffects(genotFitness = df3x,
390
-                            frequencyDependentFitness = TRUE, frequencyType = "abs")
385
+                              frequencyDependentFitness = TRUE, frequencyType = "abs")
391 386
 
392 387
     interventions = list(
393 388
         list(
394 389
             ID            = "intOverA",
395 390
             Trigger       = "T >= 1.2",
396
-            WhatHappens   = "n_A = n_A * 0.5",
391
+            WhatHappens   = "n_A = n_A * 0.2",
397 392
             Repetitions   = 0,
398 393
             Periodicity   = Inf
399 394
         ),
400 395
         list(
401 396
             ID            = "intOverB",
402 397
             Trigger       = "T >= 2.2",
403
-            WhatHappens   = "n_B = n_B * 0.5",
398
+            WhatHappens   = "n_B = n_B * 0.2",
404 399
             Repetitions   = 0,
405 400
             Periodicity   = Inf
406 401
         ),
407 402
         list(
408 403
             ID            = "intOverC",
409 404
             Trigger       = "T >= 3.2",
410
-            WhatHappens   = "n_C = n_C * 0.5",
405
+            WhatHappens   = "n_C = n_C * 0.1",
411 406
             Repetitions   = 0,
412 407
             Periodicity   = Inf
413 408
         ),
414 409
         list(
415 410
             ID            = "intOverD",
416 411
             Trigger       = "T >= 4.2",
417
-            WhatHappens   = "n_D = n_D * 0.5",
412
+            WhatHappens   = "n_D = n_D * 0.01",
418 413
             Repetitions   = 0,
419 414
             Periodicity   = Inf
420 415
         )
421 416
     )
422 417
 
423 418
     interventions = createInterventions(interventions, afd3)
424
-
419
+    seed <- round(runif(1, 1, 1e9))
420
+    set.seed(seed)
425 421
     sfd3_with_ints <- oncoSimulIndiv( afd3,
426
-                            model = "McFLD",
427
-                            onlyCancer = FALSE,
428
-                            finalTime = 6,
429
-                            mu = 1e-4,
430
-                            sampleEvery = 0.01,
431
-                            keepPhylog = FALSE,
432
-                            seed = NULL,
433
-                            errorHitMaxTries = FALSE,
434
-                            errorHitWallTime = FALSE,
435
-                            initMutant = c("A", "B", "C", "D", "E"),
436
-                            initSize = c(20000, 20000, 30, 10, 200),
437
-                            interventions = interventions)
438
-
439
-    indexes <- vector()
440
-    # We get the indexes that match the intervention times in pops.by.time
441
-    for(time in sfd3_with_ints$other$interventionTimes){
442
-        indexes <- append(indexes, which(sfd3_with_ints$pops.by.time[,1] == time))
443
-    }
444
-
445
-    # when T = 1.2 population of genotype A is intervened
446
-    line <- sfd3_with_ints$pops.by.time[indexes[1], ]
447
-    prev_line <- sfd3_with_ints$pops.by.time[indexes[1]-1, ]
448
-    if((prev_line[2] > 0) & (line[2] > 0)){
449
-        testthat::expect_gt(line[2], prev_line[2]*0.5 - 0.2*prev_line[2])
450
-        testthat::expect_lt(line[2], prev_line[2]*0.5 + 0.2*prev_line[2])
451
-    }
452
-
453
-    # when T = 2.2 population of genotype B is intervened
454
-    line <- sfd3_with_ints$pops.by.time[indexes[2], ]
455
-    prev_line <- sfd3_with_ints$pops.by.time[indexes[2]-1, ]
456
-    if((prev_line[3] > 0) & (line[3] > 0)){
457
-        testthat::expect_gt(line[3], prev_line[3]*0.5 - 0.2*prev_line[3])
458
-        testthat::expect_lt(line[3], prev_line[3]*0.5 + 0.2*prev_line[3])
459
-    }
460
-    # when T = 3.2 population of genotype C is intervened
461
-    line <- sfd3_with_ints$pops.by.time[indexes[3], ]
462
-    prev_line <- sfd3_with_ints$pops.by.time[indexes[3]-1, ]
463
-    if((prev_line[4] > 0) & (line[4] > 0)){
464
-        testthat::expect_gt(line[4], prev_line[4]*0.5 - 0.2*prev_line[4])
465
-        testthat::expect_lt(line[4], prev_line[4]*0.5 + 0.2*prev_line[4])
466
-    }
467
-    # when T = 4.2 population of genotype D is intervened
468
-    line <- sfd3_with_ints$pops.by.time[indexes[4], ]
469
-    prev_line <- sfd3_with_ints$pops.by.time[indexes[4]-1, ]
470
-    if((prev_line[5] > 0) & (line[5] > 0)){
471
-        testthat::expect_gt(line[5], prev_line[5]*0.5 - 0.2*prev_line[5])
472
-        testthat::expect_lt(line[5], prev_line[5]*0.5 + 0.2*prev_line[5])
473
-    }
422
+                                     model = "Exp",
423
+                                     onlyCancer = FALSE,
424
+                                     finalTime = 5,
425
+                                     mu = 1e-4,
426
+                                     sampleEvery = 1,
427
+                                     keepPhylog = FALSE,
428
+                                     errorHitMaxTries = FALSE,
429
+                                     errorHitWallTime = FALSE,
430
+                                     initMutant = c("A", "B", "C", "D", "E"),
431
+                                     initSize = c(20000, 20000, 3000, 10000, 200),
432
+                                     interventions = interventions)
433
+    set.seed(seed)
434
+    sfd3_without_ints <- oncoSimulIndiv( afd3,
435
+                                        model = "Exp",
436
+                                        onlyCancer = FALSE,
437
+                                        finalTime = 5,
438
+                                        mu = 1e-4,
439
+                                        sampleEvery = 1,
440
+                                        keepPhylog = FALSE,
441
+                                        errorHitMaxTries = FALSE,
442
+                                        errorHitWallTime = FALSE,
443
+                                        initMutant = c("A", "B", "C", "D", "E"),
444
+                                        initSize = c(20000, 20000, 3000, 10000, 200))
445
+    
446
+    indexes <- which(sfd3_with_ints$pops.by.time[,1] %in%
447
+                     sfd3_with_ints$other$interventionTimes)
448
+
449
+
450
+    ## Same logic as above: compare populations without and without intervention
451
+    ## But after the first difference, trajectories could differ from
452
+    ## random number divergences
453
+    expect_equal(sfd3_with_ints$pops.by.time[indexes[1], 2],
454
+                 floor(.2 * sfd3_without_ints$pops.by.time[indexes[1], 2]))
455
+    expect_lt(sfd3_with_ints$pops.by.time[indexes[2], 3],
456
+              sfd3_without_ints$pops.by.time[indexes[2], 3])
457
+    expect_lt(sfd3_with_ints$pops.by.time[indexes[3], 4],
458
+              sfd3_without_ints$pops.by.time[indexes[3], 4])
459
+    expect_lt(sfd3_with_ints$pops.by.time[indexes[4], 5],
460
+              sfd3_without_ints$pops.by.time[indexes[4], 5])
461
+
462
+    
474 463
 })
475 464
 
476 465
 
477
-
478
-test_that("12. Intervening over 4 genotypes both over specific and total population (Exp) | Trigger depends on T", {
479
-
466
+test_that("12. Intervening over 4 genotypes both over specific genotype and total population (McFL) | Trigger depends on N",{
467
+    ## Note same logic as in 11, because death rate in McFL decreases when we decrease
468
+    ## population size. So intervene at exact same time on all
469
+    ## and test exact proportions for all
480 470
     df3x <- data.frame(Genotype = c("A", "B", "C", "D", "E"),
481
-                      Fitness = c("1",
482
-                                  "1.01 + (0 * n_A)",
483
-                                  "1.1",
484
-                                  "1.09",
485
-                                  "1.07"))
471
+                       Fitness = c("1",
472
+                                   "1.01 + (0 * n_A)",
473
+                                   "1.1",
474
+                                   "1.09",
475
+                                   "1.07"))
486 476
 
487 477
     afd3 <- allFitnessEffects(genotFitness = df3x,
488
-                            frequencyDependentFitness = TRUE, frequencyType = "abs")
478
+                              frequencyDependentFitness = TRUE, frequencyType = "abs")
489 479
 
490 480
     interventions = list(
491 481
         list(
492 482
             ID            = "intOverA",
493 483
             Trigger       = "T >= 1.2",
494
-            WhatHappens   = "n_A = n_A * 0.5",
484
+            WhatHappens   = "n_A = n_A * 0.2",
495 485
             Repetitions   = 0,
496 486
             Periodicity   = Inf
497 487
         ),
498 488
         list(
499 489
             ID            = "intOverB",
500
-            Trigger       = "T >= 2.2",
501
-            WhatHappens   = "n_B = n_B * 0.5",
490
+            Trigger       = "T >= 1.2",
491
+            WhatHappens   = "n_B = n_B * 0.42",
502 492
             Repetitions   = 0,
503 493
             Periodicity   = Inf
504 494
         ),
505 495
         list(
506 496
             ID            = "intOverC",
507
-            Trigger       = "T >= 3.2",
508
-            WhatHappens   = "n_C = n_C * 0.5",
497
+            Trigger       = "T >= 1.2",
498
+            WhatHappens   = "n_C = n_C * 0.35",
509 499
             Repetitions   = 0,
510 500
             Periodicity   = Inf
511 501
         ),
512 502
         list(
513 503
             ID            = "intOverD",
514
-            Trigger       = "T >= 4.2",
515
-            WhatHappens   = "n_D = n_D * 0.5",
504
+            Trigger       = "T >= 1.2",
505
+            WhatHappens   = "n_D = n_D * 0.125",
516 506
             Repetitions   = 0,
517 507
             Periodicity   = Inf
518 508
         )
519 509
     )
520 510
 
521 511
     interventions = createInterventions(interventions, afd3)
522
-
512
+    seed <- round(runif(1, 1, 1e9))
513
+    set.seed(seed)
523 514
     sfd3_with_ints <- oncoSimulIndiv( afd3,
524
-                            model = "Exp",
525
-                            onlyCancer = FALSE,
526
-                            finalTime = 6,
527
-                            mu = 1e-4,
528
-                            sampleEvery = 0.01,
529
-                            keepPhylog = FALSE,
530
-                            seed = NULL,
531
-                            errorHitMaxTries = FALSE,
532
-                            errorHitWallTime = FALSE,
533
-                            initMutant = c("A", "B", "C", "D", "E"),
534
-                            initSize = c(20000, 20000, 30, 10, 200),
535
-                            interventions = interventions)
536
-
537
-    #plot(sfd3, show = "genotypes", type = "line")
538
-
539
-     indexes <- vector()
540
-    # We get the indexes that match the intervention times in pops.by.time
541
-    for(time in sfd3_with_ints$other$interventionTimes){
542
-        indexes <- append(indexes, which(sfd3_with_ints$pops.by.time[,1] == time))
543
-    }
544
-
545
-    # when T = 1.2 population of genotype A is intervened
546
-    line <- sfd3_with_ints$pops.by.time[indexes[1], ]
547
-    prev_line <- sfd3_with_ints$pops.by.time[indexes[1]-1, ]
548
-    if((prev_line[2] > 0) & (line[2] > 0)){
549
-        testthat::expect_gt(line[2], prev_line[2]*0.5 - 0.2*prev_line[2])
550
-        testthat::expect_lt(line[2], prev_line[2]*0.5 + 0.2*prev_line[2])
551
-    }
552
-
553
-    # when T = 2.2 population of genotype B is intervened
554
-    line <- sfd3_with_ints$pops.by.time[indexes[2], ]
555
-    prev_line <- sfd3_with_ints$pops.by.time[indexes[2]-1, ]
556
-    if((prev_line[3] > 0) & (line[3] > 0)){
557
-        testthat::expect_gt(line[3], prev_line[3]*0.5 - 0.2*prev_line[3])
558
-        testthat::expect_lt(line[3], prev_line[3]*0.5 + 0.2*prev_line[3])
559
-    }
560
-    # when T = 3.2 population of genotype C is intervened
561
-    line <- sfd3_with_ints$pops.by.time[indexes[3], ]
562
-    prev_line <- sfd3_with_ints$pops.by.time[indexes[3]-1, ]
563
-    if((prev_line[4] > 0) & (line[4] > 0)){
564
-        testthat::expect_gt(line[4], prev_line[4]*0.5 - 0.2*prev_line[4])
565
-        testthat::expect_lt(line[4], prev_line[4]*0.5 + 0.2*prev_line[4])
566
-    }
567
-    # when T = 4.2 population of genotype D is intervened
568
-    line <- sfd3_with_ints$pops.by.time[indexes[4], ]
569
-    prev_line <- sfd3_with_ints$pops.by.time[indexes[4]-1, ]
570
-    if((prev_line[5] > 0) & (line[5] > 0)){
571
-        testthat::expect_gt(line[5], prev_line[5]*0.5 - 0.2*prev_line[5])
572
-        testthat::expect_lt(line[5], prev_line[5]*0.5 + 0.2*prev_line[5])
573
-    }
515
+                                     model = "McFL",
516
+                                     onlyCancer = FALSE,
517
+                                     finalTime = 1.3,
518
+                                     mu = 1e-4,
519
+                                     sampleEvery = 0.025,
520
+                                     keepPhylog = FALSE,
521
+                                     errorHitMaxTries = FALSE,
522
+                                     errorHitWallTime = FALSE,
523
+                                     initMutant = c("A", "B", "C", "D", "E"),
524
+                                     initSize = c(20000, 20000, 30, 10, 200),
525
+                                     interventions = interventions)
526
+    set.seed(seed)
527
+    sfd3_without_ints <- oncoSimulIndiv( afd3,
528
+                                        model = "McFL",
529
+                                        onlyCancer = FALSE,
530
+                                        finalTime = 1.3,
531
+                                        mu = 1e-4,
532
+                                        sampleEvery = 0.025,
533
+                                        keepPhylog = FALSE,
534
+                                        errorHitMaxTries = FALSE,
535
+                                        errorHitWallTime = FALSE,
536
+                                        initMutant = c("A", "B", "C", "D", "E"),
537
+                                        initSize = c(20000, 20000, 30, 10, 200))
538
+    
539
+    indexes <- which(sfd3_with_ints$pops.by.time[,1] %in%
540
+                     sfd3_with_ints$other$interventionTimes)
541
+
542
+
543
+    expect_equal(sfd3_with_ints$pops.by.time[indexes, 2],
544
+                 floor(.2 * sfd3_without_ints$pops.by.time[indexes, 2]))
545
+    expect_equal(sfd3_with_ints$pops.by.time[indexes, 3],
546
+                 floor(.42 * sfd3_without_ints$pops.by.time[indexes, 3]))
547
+    expect_equal(sfd3_with_ints$pops.by.time[indexes, 4],
548
+                 floor(.35 * sfd3_without_ints$pops.by.time[indexes, 4]))
549
+    expect_equal(sfd3_with_ints$pops.by.time[indexes, 5],
550
+                 floor(.125 * sfd3_without_ints$pops.by.time[indexes, 5]))
551
+    
574 552
 })
575 553
 
576 554
 
577 555
 
578 556
 
579
-test_that("13. Intervening in the Rock-Paper-Scissors model for bacterial community by Kerr.", {
580
-    crs <- function (a, b, c){
581
-        data.frame(Genotype = c("WT", "C", "R"),
582
-        Fitness = c(paste0("1 + ", a, " * n_R/N - ", b, " * n_C/N"),
583
-        paste0("1 + ", b, " * n_/N - ", c, " * n_R/N"),
584
-        paste0("1 + ", c, " * n_C/N - ", a, " * n_/N")
585
-        ))
586
-    }
587
-
588
-    afcrs1 <- allFitnessEffects(genotFitness = crs(1, 1, 1),
589
-    frequencyDependentFitness = TRUE,
590
-    frequencyType = "abs")
591
-
592
-    lista_intervenciones = list(
593
-        list(
594
-            ID = "Bothering R strain, by reducing C",
595
-            Trigger = "n_C >= 500",
596
-            WhatHappens = "n_C = n_C * 0.1",
597
-            Periodicity = 3,
598
-            Repetitions = Inf
599
-        )
600
-    )
601
-
602
-    final_interventions = createInterventions(interventions = lista_intervenciones, genotFitness = afcrs1)
603
-
604
-    resultscrs1 <- oncoSimulIndiv(afcrs1,
605
-                                model = "McFL",
606
-                                onlyCancer = FALSE,
607
-                                finalTime = 100,
608
-                                mu = 1e-2,
609
-                                initSize = 4000,
610
-                                keepPhylog = TRUE,
611
-                                seed = NULL,
612
-                                errorHitMaxTries = FALSE,
613
-                                errorHitWallTime = FALSE,
614
-                                interventions = final_interventions)
615
-
616
-    # by reducing C, R wont spread in the population. This will mean that, with the apropiate
617
-    # periodicty in the intervention, C will never surpass WT. We try to check this in these tests.
618
-    i <- 1
619
-    while(i <= nrow(resultscrs1$pops.by.time)){
620
-        cur_nWT = resultscrs1$pops.by.time[i, 2:2]
621
-        cur_nR = resultscrs1$pops.by.time[i, 4:4]
622
-        testthat::expect_gt(cur_nWT, cur_nR)
623
-        i <- i + 1
624
-    }
625
-})
557
+## ## This is  not testing interventions, but interventions
558
+## ## plus other things in a complex to reason about model
559
+## test_that("13. Intervening in the Rock-Paper-Scissors model for bacterial community by Kerr.", {
560
+##     crs <- function (a, b, c){
561
+##         data.frame(Genotype = c("WT", "C", "R"),
562
+##         Fitness = c(paste0("1 + ", a, " * n_R/N - ", b, " * n_C/N"),
563
+##         paste0("1 + ", b, " * n_/N - ", c, " * n_R/N"),
564
+##         paste0("1 + ", c, " * n_C/N - ", a, " * n_/N")
565
+##         ))
566
+##     }
567
+
568
+##     afcrs1 <- allFitnessEffects(genotFitness = crs(1, 1, 1),
569
+##     frequencyDependentFitness = TRUE,
570
+##     frequencyType = "abs")
571
+
572
+##     lista_intervenciones = list(
573
+##         list(
574
+##             ID = "Bothering R strain, by reducing C",
575
+##             Trigger = "n_C >= 500",
576
+##             WhatHappens = "n_C = n_C * 0.1",
577
+##             Periodicity = 3,
578
+##             Repetitions = Inf
579
+##         )
580
+##     )
581
+
582
+##     final_interventions = createInterventions(interventions = lista_intervenciones, genotFitness = afcrs1)
583
+
584
+##     resultscrs1 <- oncoSimulIndiv(afcrs1,
585
+##                                 model = "McFL",
586
+##                                 onlyCancer = FALSE,
587
+##                                 finalTime = 100,
588
+##                                 mu = 1e-2,
589
+##                                 initSize = 4000,
590
+##                                 keepPhylog = FALSE,
591
+##                                 seed = NULL,
592
+##                                 errorHitMaxTries = FALSE,
593
+##                                 errorHitWallTime = FALSE,
594
+##                                 interventions = final_interventions)
595
+
596
+##     # by reducing C, R wont spread in the population. This will mean that, with the apropiate
597
+##     # periodicty in the intervention, C will never surpass WT. We try to check this in these tests.
598
+##     i <- 1
599
+##     while(i <= nrow(resultscrs1$pops.by.time)){
600
+##         cur_nWT = resultscrs1$pops.by.time[i, 2:2]
601
+##         cur_nR = resultscrs1$pops.by.time[i, 4:4]
602
+##         testthat::expect_gt(cur_nWT, cur_nR)
603
+##         i <- i + 1
604
+##     }
605
+## })
626 606
 
627 607
 test_that("14. Intervening over total population (Exp) | Trigger depends on user variables", {
628 608
     gffd3 <- data.frame(Genotype = c("WT", "A", "B"),
629
-                    Fitness = c("1",
630
-                    "1 + 0.2 * (n_B > 0)",
631
-                    ".9 + 0.4 * (n_A > 0)"
632
-                    ))
609
+                        Fitness = c("1",
610
+                                    "1 + 0.2 * (n_B > 0)",
611
+                                    ".9 + 0.4 * (n_A > 0)"
612
+                                    ))
633 613
     afd3 <- allFitnessEffects(genotFitness = gffd3,
634
-                            frequencyDependentFitness = TRUE,
635
-                            frequencyType = "abs")
614
+                              frequencyDependentFitness = TRUE,
615
+                              frequencyType = "abs")
636 616
 
637 617
     userVars <- list(
638 618
         list(Name = "user_var_1",
639
-            Value = 0
640
-        )
619
+             Value = 0
620
+             )
641 621
     )
642 622
 
643 623
     userVars <- createUserVars(userVars)
644 624
 
645 625
     rules <- list(
646 626
         list(ID = "rule_1",
647
-            Condition = "T >= 10",
648
-            Action = "user_var_1 = 1"
649
-        ),list(ID = "rule_2",
650
-            Condition = "T >= 20",
651
-            Action = "user_var_1 = 2"
652
-        ),list(ID = "rule_3",
653
-            Condition = "T >= 30",
654
-            Action = "user_var_1 = 3"
655
-        )
627
+             Condition = "T >= 10",
628
+             Action = "user_var_1 = 1"
629
+             ),list(ID = "rule_2",
630
+                    Condition = "T >= 20",
631
+                    Action = "user_var_1 = 2"
632
+                    ),list(ID = "rule_3",
633
+                           Condition = "T >= 30",
634
+                           Action = "user_var_1 = 3"
635
+                           )
656 636
     )
657 637
 
658 638
     rules <- createRules(rules, afd3)
... ...
@@ -660,112 +640,84 @@ test_that("14. Intervening over total population (Exp) | Trigger depends on user
660 640
         list(
661 641
             ID            = "intOverTotPop",
662 642
             Trigger       = "user_var_1 = 1",
663
-            WhatHappens   = "N = N * 0.8",
664
-            Repetitions   = 0,
665
-            Periodicity   = Inf
666
-        ),list(
667
-            ID            = "intOverTotPop2",
668
-            Trigger       = "user_var_1 = 2",
669
-            WhatHappens   = "N = N * 0.8",
643
+            WhatHappens   = "N = N * 0.7",
670 644
             Repetitions   = 0,
671 645
             Periodicity   = Inf
672 646
         ),list(
673
-            ID            = "intOverTotPop3",
674
-            Trigger       = "user_var_1 = 3",
675
-            WhatHappens   = "N = N * 0.8",
676
-            Repetitions   = 0,
677
-            Periodicity   = Inf
678
-        )
647
+              ID            = "intOverTotPop2",
648
+              Trigger       = "user_var_1 = 2",
649
+              WhatHappens   = "N = N * 0.2",
650
+              Repetitions   = 0,
651
+              Periodicity   = Inf
652
+          ),list(
653
+                ID            = "intOverTotPop3",
654
+                Trigger       = "user_var_1 = 3",
655
+                WhatHappens   = "N = N * 0.5",
656
+                Repetitions   = 0,
657
+                Periodicity   = Inf
658
+            )
679 659
     )
680 660
 
681 661
     interventions <- createInterventions(interventions, afd3)
682 662