tests/testthat/test-hmmbasic-riself4-8-16.R

In qtl2: Quantitative Trait Locus Mapping in Experimental Crosses

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)

})