Update test for CS_version = 0

rsetienne · May 22, 2021 · e05e038 · e05e038
1 parent 842f5b3
commit e05e038
Show file tree

Hide file tree

Showing 3 changed files with 72 additions and 16 deletions.
diff --git a/R/DAISIE_loglik_CS.R b/R/DAISIE_loglik_CS.R
@@ -835,8 +835,7 @@ DAISIE_loglik_CS <- DAISIE_loglik_all <- function(pars1,
         methode = methode,
         CS_version = CS_version,
         abstolint = abstolint,
-        reltolint = reltolint
-      )
+        reltolint = reltolint)
     }
   }
   if (pars2[4] >= 1) {

diff --git a/R/DAISIE_loglik_IW0.R b/R/DAISIE_loglik_IW0.R
@@ -5,7 +5,7 @@ DAISIE_loglik_IW_M1 <- function(
   brts,
   stac,
   missnumspec,
-  methode = "ode45",
+  methode,
   abstolint = 1E-16,
   reltolint = 1E-14,
   verbose
@@ -14,19 +14,33 @@ DAISIE_loglik_IW_M1 <- function(
     stop('This likelihood computation cannot deal with missing species.')
   }
   if(!(stac %in% c(0,2,4))) {
-    stop('This likelihood computation must have explicit colonization times.')
+    stop('This likelihood computation must have explicit colonization times or none at all.')
   }
   datalist2 <- list()
   datalist2[[1]] <- list(island_age = max(abs(brts)), not_present = as.numeric(is.null(datalist)))
-  if(!is.null(datalist)) datalist2[[2]] <- datalist
-  loglik <- DAISIE_loglik_IW(
-    pars1 = pars1,
-    pars2 = pars2,
-    datalist = datalist2,
-    methode = 'odeint::runge_kutta_fehlberg78',
-    abstolint = abstolint,
-    reltolint = reltolint,
-    verbose = verbose
-  )
+  if(!is.null(datalist)) {
+    pars2[3] <- 0
+    datalist2[[2]] <- datalist
+    loglik <- DAISIE_loglik_IW(
+      pars1 = pars1,
+      pars2 = pars2,
+      datalist = datalist2,
+      methode = methode,
+      abstolint = abstolint,
+      reltolint = reltolint,
+      verbose = verbose)
+  } else {
+    loglik <- DAISIE_loglik(
+      pars1 = pars1,
+      pars2 = pars2,
+      brts = max(abs(brts)),
+      stac = 0,
+      missnumspec = missnumspec,
+      methode = 'ode45',
+      abstolint = abstolint,
+      reltolint = reltolint,
+      verbose = verbose
+    )
+  }
   return(loglik)
 }
diff --git a/tests/testthat/test-DAISIE_loglik_CS.R b/tests/testthat/test-DAISIE_loglik_CS.R
@@ -45,7 +45,7 @@ test_that("DAISIE_loglik_CS_choice produces correct output for relaxed-rate
 
 test_that("DAISIE_loglik_CS_choice produces same output for CS_version = 0 (with M = 1) and CS_version = 1 ", {
   pars1 <- c(2.000, 2.700, 20.000, 0.009, 1.010)
-  pars2 <- c(1.0e+02, 1.1e+01, 0.0e+00, 0.0e+00, NA, 0.0e+00, 1.0e-04,
+  pars2 <- c(100, 11, 0, 0, NA, 0.0e+00, 1.0e-04,
              1.0e-05, 1.0e-07, 3.0e+03, 9.5e-01, 9.8e-01)
   brts <- c(4.0000, 3.0282, 1.3227, 0.8223, 0.4286, 0.3462, 0.2450, 0.0808,
             0.0527, 0.0327, 0.0221, 0.1180, 0.0756, 0.0525, 0.0322, 0.0118)
@@ -69,7 +69,6 @@ test_that("DAISIE_loglik_CS_choice produces same output for CS_version = 0 (with
                                     CS_version = CS_version)
 
   expect_equal(loglik0,loglik1)
-
 })
 
 test_that("DAISIE_loglik_all produces correct output for relaxed-rate model", {
@@ -104,3 +103,47 @@ test_that("DAISIE_loglik produces correct output", {
                             verbose = FALSE)
   testthat::expect_equal(output, -0.00347317077256095)
 })
+
+test_that("DAISIE_loglik_all produces same output for CS_version 0 and 1 with and without conditioning", {
+  utils::data(Galapagos_datalist)
+  Galapagos_datalist2 <- Galapagos_datalist
+  for(i in 2:9) {
+    Galapagos_datalist2[[i]]$branching_times <- c(4, 4 - 2*i*0.1,4 -2*i*0.1-0.1)
+    Galapagos_datalist2[[i]]$stac <- 2
+  }
+  Galapagos_datalist2 <- DAISIE:::add_brt_table(Galapagos_datalist2)
+  loglik_CS00 <- DAISIE::DAISIE_loglik_all(
+    pars1 = c(2.55068735, 2.68345455, 10.00000000, 0.00933207, 1.01007312),
+    pars2 = c(100, 11, 0, 0, NA),
+    datalist = Galapagos_datalist2,
+    methode = "odeint::runge_kutta_fehlberg78",
+    CS_version = 0,
+    abstolint = 1e-16,
+    reltolint = 1e-10)
+  loglik_CS10 <- DAISIE::DAISIE_loglik_all(
+    pars1 = c(2.55068735, 2.68345455, 10.00000000, 0.00933207, 1.01007312),
+    pars2 = c(100, 11, 0, 0, NA),
+    datalist = Galapagos_datalist2,
+    methode = "ode45",
+    CS_version = 1,
+    abstolint = 1e-16,
+    reltolint = 1e-10)
+  testthat::expect_equal(loglik_CS00, loglik_CS10, tol = 5E-6)
+  loglik_CS01 <- DAISIE::DAISIE_loglik_all(
+    pars1 = c(2.55068735, 2.68345455, 10.00000000, 0.00933207, 1.01007312),
+    pars2 = c(100, 11, 1, 0, NA),
+    datalist = Galapagos_datalist2,
+    methode = "odeint::runge_kutta_fehlberg78",
+    CS_version = 0,
+    abstolint = 1e-16,
+    reltolint = 1e-10)
+  loglik_CS11 <- DAISIE::DAISIE_loglik_all(
+    pars1 = c(2.55068735, 2.68345455, 10.00000000, 0.00933207, 1.01007312),
+    pars2 = c(100, 11, 1, 0, NA),
+    datalist = Galapagos_datalist2,
+    methode = "ode45",
+    CS_version = 1,
+    abstolint = 1e-16,
+    reltolint = 1e-10)
+  testthat::expect_equal(loglik_CS01, loglik_CS11, tol = 5E-6)
+})