diff --git a/R/DAISIE_format_IW.R b/R/DAISIE_format_IW.R index 2af0a72a..c26d1867 100644 --- a/R/DAISIE_format_IW.R +++ b/R/DAISIE_format_IW.R @@ -117,7 +117,7 @@ DAISIE_format_IW_trait <- function(island_replicates, return(several_islands) } -add_brt_table <- function(island, full_table = FALSE) { +add_brt_table_old <- function(island, full_table = FALSE) { island_age <- island[[1]]$island_age island_top <- island[[1]] if (length(island) == 1) { @@ -204,8 +204,7 @@ add_brt_table <- function(island, full_table = FALSE) { return(island) } -#The function below makes more sense, but DAISIE_loglik_IW needs to be adjusted to it -add_brt_table2 <- function(island, full_table = FALSE) { +add_brt_table <- function(island, full_table = TRUE) { island_age <- island[[1]]$island_age island_top <- island[[1]] if (length(island) == 1) { @@ -247,7 +246,7 @@ add_brt_table2 <- function(island, full_table = FALSE) { for (i in 1:length(btimes)) { the_stac <- island_no_stac1or5[[i]]$stac if(!is.null(island[[i]]$all_colonisations) & full_table == TRUE) { - if(length(island[[i]]$all_colonisations) > 0) print(i) + #if(length(island[[i]]$all_colonisations) > 0) print(i) for (k in 1:length(island[[i]]$all_colonisations)) { the_brts <- island[[i]]$all_colonisations[[k]]$event_times[-1] pos2 <- pos1 + length(the_brts) diff --git a/R/DAISIE_loglik_IW.R b/R/DAISIE_loglik_IW.R index 6a059c48..099be6e6 100644 --- a/R/DAISIE_loglik_IW.R +++ b/R/DAISIE_loglik_IW.R @@ -427,7 +427,12 @@ DAISIE_loglik_IW <- function( test_for_colonization <- TRUE } else { - test_for_colonization <- (max(event[which(clade == col[k + 2])]) > 1) + #test_for_colonization <- (max(event[which(clade == col[k + 2])]) > 1) + # this tests whether the original clade has diversified, but this does not need to be the case + # it could also be a case of anagenesis followed by colonization, so this test seems inappropriate + # also, one should expect that colonizations in the data are all real colonizations; a colonization + # which is not followed by speciation but by recolonization should only occur once (the last colonization) + test_for_colonization <- TRUE } if(test_for_colonization) # new colonization or recolonization after speciation { diff --git a/tests/testthat/test-DAISIE_format_IW.R b/tests/testthat/test-DAISIE_format_IW.R index d257aa3a..fcdfab0e 100644 --- a/tests/testthat/test-DAISIE_format_IW.R +++ b/tests/testthat/test-DAISIE_format_IW.R @@ -94,11 +94,11 @@ test_that("silent with non-empty island with correct output", { brts_table <- matrix(ncol = 5, nrow = 6) colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col") brts_table[1, ] <- c(1, 0, 0, NA, NA) - brts_table[2, ] <- c(0.9244818166871660, 1, 1, 1, NA) - brts_table[3, ] <- c(0.9105856673960619, 1, 2, 1, NA) - brts_table[4, ] <- c(0.5557734125062590, 2, 1, 0, NA) - brts_table[5, ] <- c(0.5288428248966160, 3, 1, 0, NA) - brts_table[6, ] <- c(0.3146835586399670, 1, 3, 1, NA) + brts_table[2, ] <- c(0.9244818166871660, 1, 1, 1, 1) + brts_table[3, ] <- c(0.9105856673960619, 1, 2, 1, 1) + brts_table[4, ] <- c(0.5557734125062590, 2, 1, 0, 1) + brts_table[5, ] <- c(0.5288428248966160, 3, 1, 0, 1) + brts_table[6, ] <- c(0.3146835586399670, 1, 3, 1, 1) expected_IW_format[[1]][[1]] <- list(island_age = 1, not_present = 7, stt_all = stt_all, @@ -255,8 +255,8 @@ test_that("silent with non-empty nonoceanic island with brts_table <- matrix(ncol = 5, nrow = 3) colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col") brts_table[1, ] <- c(1, 0, 0, NA, NA) - brts_table[2, ] <- c(1, 2, 1, 1, NA) - brts_table[3, ] <- c(1, 1, 1, 1, NA) + brts_table[2, ] <- c(1, 2, 1, 1, 1) + brts_table[3, ] <- c(1, 1, 1, 1, 1) expected_IW_format[[1]][[1]] <- list(island_age = 1, not_present = 8, stt_all = stt_all, @@ -316,8 +316,8 @@ test_that("silent with non-empty nonoceanic island with brts_table <- matrix(ncol = 5, nrow = 3) colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col") brts_table[1, ] <- c(1, 0, 0, NA, NA) - brts_table[2, ] <- c(1, 2, 1, 1, NA) - brts_table[3, ] <- c(1, 1, 1, 1, NA) + brts_table[2, ] <- c(1, 2, 1, 1, 1) + brts_table[3, ] <- c(1, 1, 1, 1, 1) expected_IW_format[[1]][[1]] <- list(island_age = 1, not_present = 8, stt_all = stt_all, @@ -385,10 +385,10 @@ test_that("add_brt_table output is correct when length(island) != 1", { brt_table <- matrix(ncol = 5, nrow = 5) colnames(brt_table) <- c("brt", "clade", "event", "endemic", "col") brt_table[1, ] <- c(1, 0, 0, NA, NA) - brt_table[2, ] <- c(0.9244818, 1, 1, 1, NA) - brt_table[3, ] <- c(0.9105857, 1, 2, 1, NA) - brt_table[4, ] <- c(0.5557734, 2, 1, 0, NA) - brt_table[5, ] <- c(0.3146836, 1, 3, 1, NA) + brt_table[2, ] <- c(0.9244818, 1, 1, 1, 1) + brt_table[3, ] <- c(0.9105857, 1, 2, 1, 1) + brt_table[4, ] <- c(0.5557734, 2, 1, 0, 1) + brt_table[5, ] <- c(0.3146836, 1, 3, 1, 1) expected_brt <- list() expected_brt[[1]] <- list(island_age = 1, not_present = 3, @@ -540,9 +540,9 @@ test_that("silent when species with two trait states with brts_table <- matrix(ncol = 5, nrow = 4) colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col") brts_table[1, ] <- c(5.00000000000000, 0, 0, NA, NA) - brts_table[2, ] <- c(3.10261367452990, 1, 1, 1, NA) - brts_table[3, ] <- c(1.50562999775257, 2, 1, 1, NA) - brts_table[4, ] <- c(1.26245655913561, 2, 2, 1, NA) + brts_table[2, ] <- c(3.10261367452990, 1, 1, 1, 1) + brts_table[3, ] <- c(1.50562999775257, 2, 1, 1, 1) + brts_table[4, ] <- c(1.26245655913561, 2, 2, 1, 1) expected_IW_format[[1]][[1]] <- list(island_age = 5, not_present = 13, stt_all = stt_all, diff --git a/tests/testthat/test-DAISIE_format_IW_full_stt.R b/tests/testthat/test-DAISIE_format_IW_full_stt.R index 6cd03312..6b5f6a31 100644 --- a/tests/testthat/test-DAISIE_format_IW_full_stt.R +++ b/tests/testthat/test-DAISIE_format_IW_full_stt.R @@ -67,7 +67,7 @@ test_that("complete stt, 1 type, no geodynamics, oceanic island, one trait state c(Time = 0.0, nI = 0.0, nA = 1.0, nC = 0.0) ) expected_brts_table <- matrix( - c(1.0, 0.24481816687165, 0, 1, 0, 1, NA, 1, NA, NA), + c(1.0, 0.24481816687165, 0, 1, 0, 1, NA, 1, NA, 1), nrow = 2 ) colnames(expected_brts_table) <- c("brt", "clade", "event", "endemic", "col") @@ -109,7 +109,7 @@ test_that("complete stt, 1 type, no geodynamics, oceanic island, one trait state ) expected_brts_table <- matrix( - c(1.0, 0.741771912202239, 0, 1, 0, 1, NA, 0, NA, NA), + c(1.0, 0.741771912202239, 0, 1, 0, 1, NA, 0, NA, 1), nrow = 2 ) colnames(expected_brts_table) <- c("brt", "clade", "event", "endemic", "col") @@ -245,7 +245,7 @@ test_that("complete stt, 1 type, geodynamics, oceanic island, one trait state ) testthat::expect_equal( formatted_IW_sim[[1]][[1]]$brts_table[5, ], - c(brt = 0.83094531417507, clade = 4, event = 1, endemic = 1, col = NA) + c(brt = 0.83094531417507, clade = 4, event = 1, endemic = 1, col = 1) ) }) diff --git a/tests/testthat/test-DAISIE_format_IW_sampled_stt.R b/tests/testthat/test-DAISIE_format_IW_sampled_stt.R index 34ff22bf..8424f9a4 100644 --- a/tests/testthat/test-DAISIE_format_IW_sampled_stt.R +++ b/tests/testthat/test-DAISIE_format_IW_sampled_stt.R @@ -59,9 +59,9 @@ test_that("sampled stt, 1 type, no geodynamics, oceanic island (same arguments a expected_brts_table <- matrix( data = c( 1.0, 0.0, 0.0, NA, NA, - 0.244818166871655, 1, 1, 1, NA, - 0.173128288990374, 1, 2, 1, NA, - 0.029668240213840, 1, 3, 1, NA + 0.244818166871655, 1, 1, 1, 1, + 0.173128288990374, 1, 2, 1, 1, + 0.029668240213840, 1, 3, 1, 1 ), ncol = 5, byrow = TRUE diff --git a/tests/testthat/test-integration_DAISIE.R b/tests/testthat/test-integration_DAISIE.R index 425744f3..d5f9ac79 100644 --- a/tests/testthat/test-integration_DAISIE.R +++ b/tests/testthat/test-integration_DAISIE.R @@ -100,6 +100,44 @@ test_that("IW and CS loglik is same when K = Inf", { testthat::expect_equal(loglik_IW, loglik_CS, tol = 5E-6) }) +test_that("IW and CS loglik is same when K = Inf", { + skip_if(Sys.getenv("CI") == "" && !(Sys.getenv("USERNAME") == "rampa"), + message = "Run only on CI") + skip_on_cran() + M <- 1000 + data(frogs_datalist, package = "DAISIE") + frogs_datalist[[1]]$not_present <- frogs_datalist[[1]]$not_present + (M - 300) + ddmodel <- 11 + initparsopt <- c(4.012298e-01,1.699521e-01,1.319595e+02,3.487955e-04) + idparsopt <- c(1,2,3,4) + parsfix <- 0 + idparsfix <- 5 + verbose <- 1 + cond <- 1 + res <- 200 + methode <- 'odeint::runge_kutta_fehlberg78' + tolint <- c(1E-16, 1E-14) + + trparsopt <- initparsopt / (1 + initparsopt) + trparsopt[which(initparsopt == Inf)] <- 1 + trparsfix <- parsfix / (1 + parsfix) + trparsfix[which(parsfix == Inf)] <- 1 + pars2 <- c(res, ddmodel, cond, verbose) + initloglik_IW <- DAISIE_loglik_IW_choosepar(trparsopt = trparsopt, + trparsfix = trparsfix, + idparsopt = idparsopt, + idparsfix = idparsfix, + M = M, + pars2 = pars2, + datalist = frogs_datalist, + methode = methode, + abstolint = tolint[1], + reltolint = tolint[2]) + + testthat::expect_equal(initloglik_IW, -215.097677998973, tol = 1E-6) +}) + + test_that("DAISIE_ML simple case works", { skip_if(Sys.getenv("CI") == "" && !(Sys.getenv("USERNAME") == "rampa"), message = "Run only on CI")