Browse code

2.99.93; fixed bugs and test errors for three-element vector

ramon diaz-uriarte (at Phelsuma) authored on 29/04/2021 23:22:43
Showing5 changed files

... ...
@@ -1,8 +1,8 @@
1 1
 Package: OncoSimulR
2 2
 Type: Package
3 3
 Title: Forward Genetic Simulation of Cancer Progression with Epistasis 
4
-Version: 2.99.92
5
-Date: 2021-04-27
4
+Version: 2.99.93
5
+Date: 2021-04-30
6 6
 Authors@R: c(
7 7
 	      person("Ramon", "Diaz-Uriarte", role = c("aut", "cre"),	
8 8
  	   		     email = "rdiaz02@gmail.com"),
... ...
@@ -232,14 +232,28 @@ rfitness <- function(g, c= 0.5,
232 232
             }
233 233
         } else { ## a length-3 scale
234 234
             if(scale[2] > scale[3]) warning("In scale, minimum fitness > wildtype")
235
+            if(scale[1] < scale[3]) warning("In scale, maximum fitness < wildtype")
235 236
             fiwt <- fi[1]
236
-            prod_above <- (scale[1] - scale[3]) / (max(fi) - fiwt)
237
-            prod_below <- (scale[3] - scale[2]) / (fiwt - min(fi))
238
-            fi_above <- which(fi >= fiwt)
239
-            fi_below <- which(fi < fiwt)
240
-            fi[fi_above] <- ((fi[fi_above] - fiwt) * prod_above) + scale[3]
241
-            fi[fi_below] <- ((fi[fi_below] - fiwt) * prod_below) + scale[3]
242
-            fi[1] <- scale[3]
237
+            new_fi <- rep(NA, length(fi))
238
+            mode(new_fi) <- "numeric"
239
+            ## If WT is min or max, there are no below or above
240
+            if(max(fi) == fiwt)
241
+                warning("WT has maximum fitness. Range will be from scale[2] to scale[3]")
242
+            if(min(fi) == fiwt)
243
+                warning("WT has minimum fitness. Range will be from scale[3] to scale[1]")
244
+            if(max(fi) > fiwt) {
245
+                prod_above <- (scale[1] - scale[3]) / (max(fi) - fiwt)
246
+                fi_above <- which(fi >= fiwt)
247
+                new_fi[fi_above]  <- ((fi[fi_above] - fiwt) * prod_above) + scale[3]
248
+            }
249
+            if(min(fi) < fiwt) {
250
+                prod_below <- (scale[3] - scale[2]) / (fiwt - min(fi))
251
+                fi_below <- which(fi < fiwt)
252
+                new_fi[fi_below] <- ((fi[fi_below] - fiwt) * prod_below) + scale[3]
253
+            }
254
+            new_fi[1] <- scale[3]
255
+            fi <- new_fi
256
+            rm(new_fi)
243 257
         }
244 258
         
245 259
         if(log) {
... ...
@@ -1,3 +1,7 @@
1
+Changes in version 2.99.93 (2021-04-30):
2
+	- Fixed bugs and improved testing of rfitness with
3
+	three-element scale vector.
4
+	
1 5
 Changes in version 2.99.92 (2021-04-27):
2 6
 	- rfitness: scale can take a three-element vector.
3 7
 	- Vignette: examples (not run) for deviations from SSWM.
... ...
@@ -115,7 +115,7 @@ additive models.}
115 115
   This option has no effect if you pass a three-element vector for
116 116
   \code{scale}. Using a three-element vector for \code{scale} is
117 117
   probably the most natural way of changing the scale and range of
118
-  fitness while setting the wildtype to value of your choice.
118
+  fitness while setting the wildtype to a value of your choice.
119 119
   
120 120
 }
121 121
 
... ...
@@ -101,96 +101,239 @@ test_that("Additive, Ising, Eggbox, Full exercising", {
101 101
 })
102 102
  
103 103
 test_that("Testing the three-element scale argument", {
104
-    for(j in 1:5) {
104
+    jj <- 0
105
+    for(j in 1:100) {
106
+    ## while(TRUE) {
107
+        ## jj <- jj + 1; cat(paste("\n ######### Doing jj = ", jj, "\n"))
105 108
         ## First, WT not equal to 1
106 109
         ## Don't use too stupid arguments
107 110
         scalev <- sort(runif(3, -2, 5), decreasing = TRUE)
108 111
         scalev <- scalev[c(1, 3, 2)]
109
-        i <- round(100 * runif(1))
112
+        
113
+        wt_min <- wt_max <- FALSE
114
+
115
+        i <- round(runif(1, 1, .Machine$integer.max - 1))
116
+        ## cat(paste("    i = "), i, "\n")
117
+
110 118
         set.seed(i)
111 119
         ra <- rfitness(7, truncate_at_0 = FALSE, wt_is_1 = "no")
112
-        set.seed(i)
113
-        rb <- rfitness(7, scale = scalev, truncate_at_0 = FALSE)
120
+        
121
+        if(max(ra[, "Fitness"]) == ra[1, "Fitness"]) {
122
+            set.seed(i)
123
+            expect_warning(rb <- rfitness(7, scale = scalev,
124
+                                          truncate_at_0 = FALSE),
125
+                           "WT has maximum fitness. Range will be from scale[2] to scale[3]",
126
+                           fixed = TRUE)
127
+            wt_max <- TRUE
128
+            wt_min <- FALSE
129
+        } else if(min(ra[, "Fitness"]) == ra[1, "Fitness"]) {
130
+            set.seed(i)
131
+            expect_warning(rb <- rfitness(7, scale = scalev,
132
+                                          truncate_at_0 = FALSE),
133
+                           "WT has minimum fitness. Range will be from scale[3] to scale[1]",
134
+                           fixed = TRUE)
135
+            wt_min <- TRUE
136
+            wt_max <- FALSE
137
+        } else {
138
+            set.seed(i)
139
+            rb <- rfitness(7, scale = scalev,
140
+                       truncate_at_0 = FALSE)
141
+        }
142
+        
143
+        rap <- ra[-1, "Fitness"][ra[-1, "Fitness"] > ra[1, "Fitness"]]
144
+        ran <- ra[-1, "Fitness"][ra[-1, "Fitness"] < ra[1, "Fitness"]]
145
+        rbp <- rb[-1, "Fitness"][rb[-1, "Fitness"] > scalev[3]]
146
+        rbn <- rb[-1, "Fitness"][rb[-1, "Fitness"] < scalev[3]]
147
+
148
+        ## Numerical issues with equality comparisons
149
+        rap <- c(ra[1, "Fitness"], rap)
150
+        rbp <- c(scalev[3], rbp)
151
+        ran <- c(ra[1, "Fitness"], ran)
152
+        rbn <- c(scalev[3], rbn)        
153
+
154
+        expect_true(abs(length(rap) - length(rbp)) < 2)
155
+        expect_true(abs(length(ran) - length(rbn)) < 2)
114 156
 
115
-        rap <- ra[, "Fitness"][ra[, "Fitness"] >= ra[1, "Fitness"]]
116
-        ran <- ra[, "Fitness"][ra[, "Fitness"] < ra[1, "Fitness"]]
117
-        rbp <- rb[, "Fitness"][rb[, "Fitness"] >= scalev[3]]
118
-        rbn <- rb[, "Fitness"][rb[, "Fitness"] < scalev[3]]
157
+        if((length(rap) > 1) && (length(rap) == length(rbp)))
158
+            expect_equivalent(cor(rap, rbp), 1)
159
+        if((length(ran) > 1) && (length(ran) == length(rbn)))
160
+            expect_equivalent(cor(ran, rbn), 1)
119 161
 
120
-        expect_equivalent(cor(rap, rbp), 1)
121
-        expect_equivalent(cor(ran, rbn), 1)
162
+        expect_true( (length(ran) + length(rap)) == (2^7 + 1) )
163
+        expect_true( (length(rbn) + length(rbp)) == (2^7 + 1) )
164
+                    
165
+        
166
+        if(!wt_max) {
167
+            expect_equivalent(max(rb[, "Fitness"]), scalev[1])
168
+        } else {
169
+            expect_equivalent(max(rb[, "Fitness"]), scalev[3])
170
+        }
171
+        if(!wt_min) {
172
+            expect_equivalent(min(rb[, "Fitness"]), scalev[2])
173
+        } else {
174
+            expect_equivalent(min(rb[, "Fitness"]), scalev[3])            
175
+        }
122 176
 
123
-        expect_equivalent(max(rb[, "Fitness"]), scalev[1])
124
-        expect_equivalent(min(rb[, "Fitness"]), scalev[2])
125 177
         expect_equivalent(rb[1, "Fitness"], scalev[3])
126 178
 
127 179
         ## X, Y, Z: X: original, Y: after scaling
128 180
         ## Just dealing with the positive
129 181
         ## Do not assume WT to be 1
130
-
131 182
         ## Write as original * slope + intercept
132 183
         ## so you can compare against
133 184
         ## lm(rbp ~ rap)
185
+
134 186
         xM <- max(ra[, "Fitness"])
135 187
         xW <- ra[1, "Fitness"]
136
-        y <- rap * (scalev[1] - scalev[3])/(xM - xW) +
137
-            (scalev[3] - xW * ((scalev[1] - scalev[3])/(xM - xW)))
188
+        if(!wt_max) {    
189
+            y <- rap * (scalev[1] - scalev[3])/(xM - xW) +
190
+                (scalev[3] - xW * ((scalev[1] - scalev[3])/(xM - xW)))
191
+            expect_equivalent(rbp, y)
192
+            rm(y)
193
+        } else {
194
+            expect_equivalent(scalev[3], rbp)
195
+        }
138 196
 
139
-        expect_equivalent(rbp, y)
140
-        
197
+        xm <- min(ra[, "Fitness"])
198
+        if(!wt_min) {
199
+            y <- ran * (scalev[3] - scalev[2])/(xW - xm) +
200
+                (scalev[3] - xW * ((scalev[3] - scalev[2])/(xW - xm)))
201
+            expect_equivalent(rbn, y)
202
+            rm(y)
203
+        } else {
204
+            expect_equivalent(scalev[3], rbn)
205
+        }
206
+
207
+        ## log transformation
141 208
         set.seed(i)
142
-        rc <- rfitness(7, scale = exp(scalev), truncate_at_0 = FALSE, log = TRUE)
209
+        suppressWarnings(
210
+            rc <- rfitness(7, scale = exp(scalev), truncate_at_0 = FALSE,
211
+                           log = TRUE))
212
+        rcp <- rc[-1, "Fitness"][rc[-1, "Fitness"] > scalev[3]]
213
+        rcn <- rc[-1, "Fitness"][rc[-1, "Fitness"] < scalev[3]]
214
+        
215
+        ## Numerical issues with equality comparisons
216
+        rcp <- c(scalev[3], rcp)
217
+        rcn <- c(scalev[3], rcn)
218
+        
219
+        expect_true(abs(length(rap) - length(rcp)) < 2)
220
+        expect_true(abs(length(ran) - length(rcn)) < 2)
221
+        
222
+        if((length(rap) > 1) && (length(rap) == length(rcp)))
223
+            expect_equivalent(cor(rap, exp(rcp)), 1)
224
+        if((length(ran) > 1) && (length(ran) == length(rcn)))
225
+            expect_equivalent(cor(ran, exp(rcn)), 1)
143 226
         
144
-        rcp <- rc[, "Fitness"][rc[, "Fitness"] >= scalev[3]]
145
-        rcn <- rc[, "Fitness"][rc[, "Fitness"] < scalev[3]]
146
-
147
-        expect_equivalent(cor(rap, exp(rcp)), 1)
148
-        expect_equivalent(cor(ran, exp(rcn)), 1)
149
-
150 227
         ## X, U, Z: X: original, U: after scaling; Z: after log
151 228
         ## Just dealing with the positive
152 229
         ## Do not assume WT to be 1
153 230
         ## can compare against lm(exp(rcp) ~ rap)
231
+        ## if(!isTRUE(all.equal(xM, xW, check.attributes = FALSE))) {
232
+        if(!wt_max) {
233
+            u <- rap * (exp(scalev[1]) - exp(scalev[3]))/(xM - xW) +
234
+                (exp(scalev[3]) - xW * ((exp(scalev[1]) - exp(scalev[3]))/(xM - xW)))
235
+            expect_equivalent(u, exp(rcp))
236
+        } else {
237
+            expect_equivalent(scalev[3], rcp)
238
+        }
154 239
 
155
-        u <- rap * (exp(scalev[1]) - exp(scalev[3]))/(xM - xW) +
156
-            (exp(scalev[3]) - xW * ((exp(scalev[1]) - exp(scalev[3]))/(xM - xW)))
157
-        expect_equivalent(u, exp(rcp))
158
-
240
+        ## if(!isTRUE(all.equal(xm, xW, check.attributes = FALSE))) {
241
+        if(!wt_min) {
242
+            v <- ran * (exp(scalev[3]) - exp(scalev[2]))/(xW - xm) +
243
+                (exp(scalev[3]) - xW * ((exp(scalev[3]) - exp(scalev[2]))/(xW - xm)))
244
+            expect_equivalent(exp(rcn), v)
245
+        } else {
246
+            expect_equivalent(scalev[3], rcn)
247
+        }
248
+        
159 249
         ## This shows, by the way, that the log-transformed ain't linear function of
160 250
         ## log of original
161
-        z <- log( (rap - xW) * ((exp(scalev[1]) - exp(scalev[3]))/(xM - xW)) + exp(scalev[3]))
162
-        expect_equivalent(z, rcp)
163
-
251
+        if(!isTRUE(all.equal(xM, xW, check.attributes = FALSE))) {
252
+            z <- log( (rap - xW) * ((exp(scalev[1]) - exp(scalev[3]))/(xM - xW)) + exp(scalev[3]))
253
+            expect_equivalent(z, rcp)
254
+        }
255
+        
164 256
         ############################################
165 257
         ## Repeat setting WT = 1
166 258
         set.seed(i)
167 259
         rax <- rfitness(7, truncate_at_0 = FALSE, wt_is_1 = "subtract")
168
-        rapx <- rax[, "Fitness"][rax[, "Fitness"] >= 1]
169
-        ranx <- rax[, "Fitness"][rax[, "Fitness"] < 1]
260
+        rapx <- rax[-1, "Fitness"][rax[-1, "Fitness"] > 1]
261
+        ranx <- rax[-1, "Fitness"][rax[-1, "Fitness"] < 1]
262
+        rapx <- c(rax[1, "Fitness"], rapx)
263
+        ranx <- c(rax[1, "Fitness"], ranx)
264
+        
265
+        expect_true(abs(length(rapx) - length(rbp)) < 2)
266
+        expect_true(abs(length(ranx) - length(rbn)) < 2)
267
+        
268
+        if((length(rapx) > 1) && (length(rapx) == length(rbp)))
269
+            expect_equivalent(cor(rapx, rbp), 1)
270
+        if((length(ranx) > 1) && (length(ranx) == length(rbn)))
271
+            expect_equivalent(cor(ranx, rbn), 1)
170 272
 
171
-        expect_equivalent(cor(rapx, rbp), 1)
172
-        expect_equivalent(cor(ranx, rbn), 1)
273
+        if(length(rapx) < 2) expect_true(length(ranx) > 2)
173 274
         
174 275
         xMx <- max(rax[, "Fitness"])
175 276
         xWx <- rax[1, "Fitness"]
176
-        yx <- rapx * (scalev[1] - scalev[3])/(xMx - xWx) +
177
-            (scalev[3] - xWx * ((scalev[1] - scalev[3])/(xMx - xWx)))
178
-
179
-        expect_equivalent(rbp, yx)
277
+        
278
+        if(!isTRUE(all.equal(xMx, xWx, check.attributes = FALSE))) {
279
+            yx <- rapx * (scalev[1] - scalev[3])/(xMx - xWx) +
280
+                (scalev[3] - xWx * ((scalev[1] - scalev[3])/(xMx - xWx)))
281
+            expect_equivalent(rbp, yx)
282
+            rm(yx)
283
+        } else {
284
+            expect_equivalent(scalev[3], rbp)
285
+        }
180 286
 
287
+        xmx <- min(rax[, "Fitness"])
288
+        if(!isTRUE(all.equal(xmx, xWx, check.attributes = FALSE))) {
289
+            yx <- ran * (scalev[3] - scalev[2])/(xWx - xmx) +
290
+                (scalev[3] - xW * ((scalev[3] - scalev[2])/(xWx - xmx)))
291
+            expect_equivalent(rbn, yx)
292
+            rm(yx)
293
+        } else {
294
+            expect_equivalent(scalev[3], rbn)
295
+        }
296
+        
297
+        expect_true(abs(length(rapx) - length(rcp)) < 2)
298
+        expect_true(abs(length(ranx) - length(rcn)) < 2)
299
+        
181 300
         ## log
182
-        expect_equivalent(cor(rapx, exp(rcp)), 1)
183
-        expect_equivalent(cor(ranx, exp(rcn)), 1)
184
-
185
-        ux <- rapx * (exp(scalev[1]) - exp(scalev[3]))/(xMx - xWx) +
186
-            (exp(scalev[3]) - xWx * ((exp(scalev[1]) - exp(scalev[3]))/(xMx - xWx)))
187
-        expect_equivalent(ux, exp(rcp))
301
+        if((length(rapx) > 1) && (length(rapx) == length(rcp)))
302
+            expect_equivalent(cor(rapx, exp(rcp)), 1)
303
+        if((length(ranx) > 1) && (length(ranx) == length(rcn)))
304
+            expect_equivalent(cor(ranx, exp(rcn)), 1)
305
+        if(!isTRUE(all.equal(xMx, xWx, check.attributes = FALSE))) {
306
+            ux <- rapx * (exp(scalev[1]) - exp(scalev[3]))/(xMx - xWx) +
307
+                (exp(scalev[3]) - xWx * ((exp(scalev[1]) - exp(scalev[3]))/(xMx - xWx)))
308
+            expect_equivalent(ux, exp(rcp))
309
+        } else {
310
+            expect_equivalent(scalev[3], rcp)
311
+        }
188 312
 
313
+         if(!isTRUE(all.equal(xm, xW, check.attributes = FALSE))) {
314
+            vx <- ran * (exp(scalev[3]) - exp(scalev[2]))/(xW - xm) +
315
+                (exp(scalev[3]) - xW * ((exp(scalev[3]) - exp(scalev[2]))/(xW - xm)))
316
+            expect_equivalent(exp(rcn), vx)
317
+        } else {
318
+            expect_equivalent(scalev[3], rcn)
319
+        }
320
+        
189 321
         ## This shows, by the way, that the log-transformed ain't linear function of
190 322
         ## log of original
191
-        zx <- log( (rapx - xWx) * ((exp(scalev[1]) - exp(scalev[3]))/(xMx - xWx)) + exp(scalev[3]))
192
-        expect_equivalent(zx, rcp)
323
+        if(!isTRUE(all.equal(xMx, xWx, check.attributes = FALSE))) {
324
+            zx <- log( (rapx - xWx) *
325
+                       ((exp(scalev[1]) - exp(scalev[3]))/(xMx - xWx)) +
326
+                       exp(scalev[3]))
327
+            expect_equivalent(zx, rcp)
328
+        } else {
329
+            expect_equivalent(scalev[3], rcp)
330
+        }
331
+        
332
+        suppressWarnings(
333
+            rm(scalev, ra, rb, rap, ran, rbp, rbn, rc, rcp, rcn, u, v, y, z,
334
+           rapx, ranx, rax, yx, ux, zx, vx, xmx, xWx, xMx, xW, xM, xm))
193 335
     }
336
+
194 337
 })
195 338
 
196 339
 cat(paste("\n Ending exercise-rfitness at", date(), "\n"))