Add a couple of demografr-induced updates and fixes

R/interface.R

@@ -55,6 +55,8 @@ population <- function(name, time, N, parent = NULL, map = FALSE,
   if (!is.character(name) || length(name) != 1)
     stop("A population name must be a character scalar value", call. = FALSE)
 
+  time <- as.integer(time)
+
   # if this population splits from a parental population, check that the parent
   # really exists and that the split time make sense given the time of appearance
   # of the parental population
@@ -71,7 +73,6 @@ population <- function(name, time, N, parent = NULL, map = FALSE,
   if (!is.null(parent) && is.logical(map) && map == FALSE)
     map <- attr(parent, "map")
 
-  time <- as.integer(time)
   if (time < 1) stop("Split time must be a non-negative integer number", call. = FALSE)
   N <- as.integer(N)
 
@@ -573,30 +574,22 @@ gene_flow <- function(from, to, rate, start, end, overlap = TRUE) {
     stop("Both 'from' and 'to' arguments must be slendr population objects",
          call. = FALSE)
 
+  # TODO: this needs some serious restructuring -- this function originated when slendr
+  # could only do spatial simulations and some consistency checks relied on spatial
+  # attributes; as a result, some checks are overlapping and/or redundant
+
   if ((has_map(from) && !has_map(to)) || (!has_map(from) && has_map(to)))
     stop("Both or neither populations must be spatial", call. = FALSE)
 
-  # make sure both participating populations are present at the start of the
-  # gene flow event (`check_present_time()` is reused from the sampling functionality)
-  # check_present_time(start, from, offset = 0)
-  # check_present_time(end, from, offset = 0)
-  # check_present_time(start, to, offset = 0)
-  # check_present_time(end, to, offset = 0)
-
-  # make sure the population is not removed during the the admixture period
-  # (note that this will only be relevant for SLiM simulations at this moment)
-  check_removal_time(start, from)
-  check_removal_time(end, from)
-  check_removal_time(start, to)
-  check_removal_time(end, to)
+  # make sure that gene flow has a sensible value between 0 and 1
+  if (rate < 0 || rate > 1)
+    stop("Gene-flow rate must be a numeric value between 0 and 1", call. = FALSE)
 
   from_name <- unique(from$pop)
   to_name <- unique(to$pop)
 
-  if (start == end)
-    stop(sprintf("No time allowed for the %s -> %s gene flow at time %s to happen",
-                 from_name, to_name, start), call. = FALSE)
-
+  # make sure that the start->end time direction of the given gene flow corresponds to
+  # the time direction implied by the split time of the populations involved
   if (from$time[1] <= start & from$time[1] <= end &
       to$time[1] <= start & to$time[1] <= end)
     comp <- `<=`
@@ -608,6 +601,31 @@ gene_flow <- function(from, to, rate, start, end, overlap = TRUE) {
 time (gene flow %s -> %s in the time window %s-%s)",
                  from_name, to_name, start, end), call. = FALSE)
 
+  if (start == end)
+    stop(sprintf("Start and end time for the %s -> %s gene flow is the same (%s)",
+                 from_name, to_name, start), call. = FALSE)
+
+  # make sure both participating populations are present at the start of the
+  # gene flow event (`check_present_time()` is reused from the sampling functionality)
+  direction <- if (start < end) "forward" else "backward"
+  tryCatch(
+    {
+      check_present_time(start, from, offset = 0, direction = direction, allow_same = TRUE)
+      check_present_time(end, from, offset = 0, direction = direction, allow_same = TRUE)
+      check_present_time(start, to, offset = 0, direction = direction, allow_same = TRUE)
+      check_present_time(end, to, offset = 0, direction = direction, allow_same = TRUE)
+
+      check_removal_time(start, from)
+      check_removal_time(end, from)
+      check_removal_time(start, to)
+      check_removal_time(end, to)
+    },
+    error = function(e) {
+      stop(sprintf("Both %s and %s must be present within the gene-flow window %s-%s",
+                   from_name, to_name, start, end), call. = FALSE)
+    }
+  )
+
   if (has_map(from) && has_map(to)) {
     # get the last specified spatial maps before the geneflow time
     region_from <- intersect_features(from[comp(from$time, start), ] %>% .[nrow(.), ])

R/utils.R

@@ -395,14 +395,14 @@ check_event_time <- function(time, pop) {
 # (used exclusively in the schedule_sampling() function to avoid situations when
 # a user would sample from a population in the same generation that it
 # would be created)
-check_present_time <- function(time, pop, offset, direction = NULL) {
+check_present_time <- function(time, pop, offset, direction = NULL, allow_same = FALSE) {
   if (is.null(direction))
     direction <- time_direction(pop)
   split_time <- get_lineage_splits(pop)[1]
 
-  if (time == split_time |
-      (direction == "backward" & time > split_time - offset) |
-      (direction == "forward" & time < split_time + offset))
+  if ((!allow_same && time == split_time) ||
+      (direction == "backward" && time > split_time - offset) ||
+      (direction == "forward" && time < split_time + offset))
     stop("Population ", pop$pop[1], " is not present at a time ", time, call. = FALSE)
 }
 

tests/testthat/test-geneflow.R

@@ -36,3 +36,29 @@ test_that("populations must be present for them to mix (non-spatial models)", {
   expect_error(gene_flow(from = pop1, to = pop2, start = 1000, end = 0, rate = 0.1),
                "Specified times are not consistent with the assumed direction", fixed = TRUE)
 })
+
+test_that("gene-flow rate must be a value between 0 and 1", {
+  pop1 <- population("pop1", N = 100, time = 100)
+  pop2 <- population("pop2", N = 100, time = 100)
+
+  error_msg <- "Gene-flow rate must be a numeric value between 0 and 1"
+  expect_error(gene_flow(from = pop1, to = pop2, start = 1000, end = 0, rate = -0.1), error_msg)
+  expect_error(gene_flow(from = pop1, to = pop2, start = 1000, end = 0, rate = 25), error_msg)
+})
+
+test_that("populations must be already created for a gene flow to happen (forward model)", {
+  pop1 <- population("pop1", N = 100, time = 100)
+  pop2 <- population("pop2", N = 100, time = 200, parent = pop1)
+
+  expect_error(
+    gene_flow(from = pop1, to = pop2, start = 50, end = 80, rate = 0.1),
+    "^Both .* and .* must be present within the gene-flow window .*-.*$"
+  )
+  expect_error(
+    gene_flow(from = pop1, to = pop2, start = 50, end = 130, rate = 0.1),
+    "Specified times are not consistent with the assumed direction"
+  )
+
+  expect_silent(gene_flow(from = pop1, to = pop2, start = 200, end = 230, rate = 0.1))
+  expect_silent(gene_flow(from = pop1, to = pop2, start = 210, end = 230, rate = 0.1))
+})