From b2c44a1e949d043ede3aaddb9a33b85a3d355186 Mon Sep 17 00:00:00 2001 From: Martin Petr Date: Thu, 30 Mar 2023 16:01:50 +0200 Subject: [PATCH] Add a couple of demografr-induced updates and fixes --- R/interface.R | 54 ++++++++++++++++++++++------------ R/utils.R | 8 ++--- tests/testthat/test-geneflow.R | 26 ++++++++++++++++ 3 files changed, 66 insertions(+), 22 deletions(-) diff --git a/R/interface.R b/R/interface.R index 54ee6f43b..a3e137130 100644 --- a/R/interface.R +++ b/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(.), ]) diff --git a/R/utils.R b/R/utils.R index 9cea275db..8d97120b1 100644 --- a/R/utils.R +++ b/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) } diff --git a/tests/testthat/test-geneflow.R b/tests/testthat/test-geneflow.R index 501e9053b..96f9c01c6 100644 --- a/tests/testthat/test-geneflow.R +++ b/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)) +})