Nothing
context("basic HMM functions in 4-, 8-, and 16-way RIL by selfing")
test_that("riself4-8-16 n_gen, n_alleles work", {
expect_equal(nalleles("riself4"), 4)
expect_equal(nalleles("riself8"), 8)
expect_equal(nalleles("riself16"), 16)
expect_equal(test_ngen("riself4", FALSE), 4)
expect_equal(test_ngen("riself8", FALSE), 8)
expect_equal(test_ngen("riself16", FALSE), 16)
})
test_that("riself4-8-16 possible_gen work", {
expect_equal(test_possible_gen("riself4", FALSE, FALSE, 1:4), 1:4)
expect_equal(test_possible_gen("riself8", FALSE, FALSE, 1:8), 1:8)
expect_equal(test_possible_gen("riself16", FALSE, FALSE, 1:16), 1:16)
})
# FIX_ME: test check_geno
test_that("riself4-8-16 init work", {
expect_equal( sapply(1:4, function(i) test_init("riself4", i, FALSE, FALSE, 1:4)), log(rep(1/4, 4)) )
expect_equal( sapply(1:8, function(i) test_init("riself8", i, FALSE, FALSE, 1:8)), log(rep(1/8, 8)) )
expect_equal( sapply(1:16, function(i) test_init("riself16", i, FALSE, FALSE, 1:16)), log(rep(1/16, 16)) )
})
# FIX_ME: test emit
test_that("riself4 step works", {
for(rf in c(0.01, 0.1, 0.45)) {
expected <- matrix(rf/(1+2*rf), ncol=4, nrow=4)
diag(expected) <- (1-rf)/(1+2*rf)
result <- matrix(ncol=4, nrow=4)
for(i in 1:4) {
for(j in 1:4) {
result[i,j] <- test_step("riself4", i, j, rf, FALSE, FALSE, 1:4)
}
}
expect_equal(result, log(expected))
}
# test_stepmatrix: output has rows that sum to 1
obs <- vapply(test_stepmatrix("riself4", c(0.01, 0.1, 0.45), FALSE, FALSE, 1:4),
function(a) rowSums(exp(a)), rep(1,4))
expect_equal(obs, matrix(1, ncol=3, nrow=4))
# test_stepmatrix: output has rows that sum to 1
obs <- vapply(test_stepmatrix("riself4", c(0.01, 0.1, 0.45), FALSE, FALSE, c(3,1,2,4)),
function(a) rowSums(exp(a)), rep(1,4))
expect_equal(obs, matrix(1, ncol=3, nrow=4))
})
test_that("riself8 step works", {
for(rf in c(0.01, 0.1, 0.45)) {
expected <- matrix(rf/2/(1+2*rf), ncol=8, nrow=8)
diag(expected) <- (1-rf)^2/(1+2*rf)
block <- list(c(1,2), c(3,4), c(5,6), c(7,8))
for(i in seq_along(block)) {
b <- block[[i]]
expected[b[1],b[2]] <- expected[b[2],b[1]] <- rf*(1-rf)/(1+2*rf)
}
result <- matrix(ncol=8, nrow=8)
for(i in 1:8) {
for(j in 1:8) {
result[i,j] <- test_step("riself8", i, j, rf, FALSE, FALSE, 1:8)
}
}
expect_equal(result, log(expected))
}
# test with different order of founders
forder <- c(3, 8, 5, 2, 1, 7, 6, 4)
forder_rev <- c(5,4,1,8,3,7,6,2)
for(rf in c(0.01, 0.1, 0.45)) {
expected <- matrix(rf/2/(1+2*rf), ncol=8, nrow=8)
diag(expected) <- (1-rf)^2/(1+2*rf)
block <- list(c(1,2), c(3,4), c(5,6), c(7,8))
for(i in seq_along(block)) {
b <- block[[i]]
expected[b[1],b[2]] <- expected[b[2],b[1]] <- rf*(1-rf)/(1+2*rf)
}
expected <- expected[forder_rev, forder_rev]
result <- matrix(ncol=8, nrow=8)
for(i in 1:8) {
for(j in 1:8) {
result[i,j] <- test_step("riself8", i, j, rf, FALSE, FALSE, forder)
}
}
expect_equal(result, log(expected))
}
# test_stepmatrix: output has rows that sum to 1
obs <- vapply(test_stepmatrix("riself8", c(0.01, 0.1, 0.45), FALSE, FALSE, 1:8),
function(a) rowSums(exp(a)), rep(1,8))
expect_equal(obs, matrix(1, ncol=3, nrow=8))
obs <- vapply(test_stepmatrix("riself8", c(0.01, 0.1, 0.45), FALSE, FALSE, forder),
function(a) rowSums(exp(a)), rep(1,8))
expect_equal(obs, matrix(1, ncol=3, nrow=8))
})
test_that("riself16 step works", {
for(rf in c(0.01, 0.1, 0.45)) {
expected <- matrix(rf/4/(1+2*rf), ncol=16, nrow=16)
expected[1:4,1:4] <- expected[5:8,5:8] <-
expected[9:12,9:12] <- expected[13:16,13:16] <- rf*(1-rf)/2/(1+2*rf)
diag(expected) <- (1-rf)^3/(1+2*rf)
block <- list(c(1,2), c(3,4), c(5,6), c(7,8),c(9,10),c(11,12),c(13,14),c(15,16))
for(i in seq_along(block)) {
b <- block[[i]]
expected[b[1],b[2]] <- expected[b[2],b[1]] <- rf*(1-rf)^2/(1+2*rf)
}
result <- matrix(ncol=16, nrow=16)
for(i in 1:16) {
for(j in 1:16) {
result[i,j] <- test_step("riself16", i, j, rf, FALSE, FALSE, 1:16)
}
}
expect_equal(result, log(expected))
}
# test with founders in a different order
forder <- c(12, 5, 15, 13, 14, 1, 2, 3, 9, 10, 6, 11, 8, 4, 16, 7)
forder_rev <- c(6,7,8,14,2,11,16,13,9,10,12,1,4,5,3,15)
for(rf in c(0.01, 0.1, 0.45)) {
expected <- matrix(rf/4/(1+2*rf), ncol=16, nrow=16)
expected[1:4,1:4] <- expected[5:8,5:8] <-
expected[9:12,9:12] <- expected[13:16,13:16] <- rf*(1-rf)/2/(1+2*rf)
diag(expected) <- (1-rf)^3/(1+2*rf)
block <- list(c(1,2), c(3,4), c(5,6), c(7,8),c(9,10),c(11,12),c(13,14),c(15,16))
for(i in seq_along(block)) {
b <- block[[i]]
expected[b[1],b[2]] <- expected[b[2],b[1]] <- rf*(1-rf)^2/(1+2*rf)
}
expected <- expected[forder_rev, forder_rev]
result <- matrix(ncol=16, nrow=16)
for(i in 1:16) {
for(j in 1:16) {
result[i,j] <- test_step("riself16", i, j, rf, FALSE, FALSE, forder)
}
}
expect_equal(result, log(expected))
}
# test_stepmatrix: output has rows that sum to 1
obs <- vapply(test_stepmatrix("riself16", c(0.01, 0.1, 0.45), FALSE, FALSE, 1:16),
function(a) rowSums(exp(a)), rep(1,16))
expect_equal(obs, matrix(1, ncol=3, nrow=16))
# test_stepmatrix: output has rows that sum to 1
obs <- vapply(test_stepmatrix("riself16", c(0.01, 0.1, 0.45), FALSE, FALSE, forder),
function(a) rowSums(exp(a)), rep(1,16))
expect_equal(obs, matrix(1, ncol=3, nrow=16))
})
test_that("riself4-8-16 geno_names work", {
expect_equal( geno_names("riself4", LETTERS[5:8], FALSE), paste0(LETTERS[5:8], LETTERS[5:8]) )
expect_equal( geno_names("riself8", LETTERS[2:9], FALSE), paste0(LETTERS[2:9], LETTERS[2:9]) )
expect_equal( geno_names("riself16", LETTERS[11:26], FALSE), paste0(LETTERS[11:26], LETTERS[11:26]) )
})
test_that("riself4-8-16 nrec work", {
x <- matrix(ncol=16, nrow=16)
x <- matrix(as.numeric(col(x) != row(x)), ncol=16)
res4 <- matrix(ncol=4, nrow=4)
res8 <- matrix(ncol=8, nrow=8)
res16 <- matrix(ncol=16, nrow=16)
for(i in 1:16) {
for(j in 1:16) {
if(i<5 && j<5) res4[i,j] <- test_nrec("riself4", i, j, FALSE, FALSE, 1:4)
if(i<9 && j<9) res8[i,j] <- test_nrec("riself8", i, j, FALSE, FALSE, 1:8)
res16[i,j] <- test_nrec("riself16", i, j, FALSE, FALSE, 1:16)
}
}
expect_equal( res4, x[1:4,1:4])
expect_equal( res8, x[1:8,1:8])
expect_equal( res16, x)
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.