Avoid silly clashes with populations named 'ancestor'

R/compilation.R

@@ -603,7 +603,7 @@ msprime <- function(model, sequence_length, recombination_rate, samples = NULL,
   if (!is.numeric(recombination_rate) | recombination_rate < 0)
     stop("Recombination rate must be a numeric value", call. = FALSE)
 
-  if (sum(model$splits$parent == "ancestor") > 1)
+  if (sum(model$splits$parent == "__pop_is_ancestor") > 1)
     stop("Multiple ancestral populations without a common ancestor would lead to\n",
          "an infinitely deep history without coalescence. Please make sure that all\n",
          "populations trace their ancestry to a single ancestral population.\n",
@@ -831,7 +831,7 @@ write_script <- function(script_target, script_source,
 compile_splits <- function(populations, generation_time, direction, end_time) {
   split_table <- lapply(populations, function(p) {
     parent <- attr(p, "parent")
-    if (is.character(parent) && parent == "ancestor") {
+    if (is.character(parent) && parent == "__pop_is_ancestor") {
       parent_name <- parent
     } else {
       parent_name <- unique(attr(p, "parent")$pop)
@@ -866,7 +866,7 @@ compile_splits <- function(populations, generation_time, direction, end_time) {
   split_table$parent_id <- lapply(
     split_table$parent,
     function(i) {
-      if (i == "ancestor") -1
+      if (i == "__pop_is_ancestor") -1
       else split_table[split_table$pop == i, ]$pop_id
     }
   ) %>% unlist()
@@ -875,7 +875,7 @@ compile_splits <- function(populations, generation_time, direction, end_time) {
   # simulation in generation 1 regardless of which "time of appearance" was
   # specified (to include it in the burnin)
   split_table %>%
-    dplyr::mutate(tsplit_gen = ifelse(parent == "ancestor", 1, tsplit_gen))
+    dplyr::mutate(tsplit_gen = ifelse(parent == "__pop_is_ancestor", 1, tsplit_gen))
 }
 
 # Process vectorized population boundaries into a table with
@@ -917,7 +917,7 @@ compile_maps <- function(populations, split_table, resolution, generation_time,
 
   # maps of ancestral populations have to be set in the first generation,
   # regardless of the specified split time
-  ancestral_pops <- split_table[split_table$parent == "ancestor", ]$pop
+  ancestral_pops <- split_table[split_table$parent == "__pop_is_ancestor", ]$pop
   ancestral_maps <- purrr::map(ancestral_pops, ~ which(map_table$pop == .x)) %>%
     purrr::map_int(~ .x[1])
   map_table[ancestral_maps, ]$time_gen <- 1
@@ -1031,7 +1031,7 @@ compile_dispersals <- function(populations, generation_time, direction,
 
   # dispersals of ancestral populations have to be set in the first generation,
   # regardless of the specified split time
-  ancestral_pops <- split_table[split_table$parent == "ancestor", ]$pop
+  ancestral_pops <- split_table[split_table$parent == "__pop_is_ancestor", ]$pop
   indices <- purrr::map(ancestral_pops, ~ which(dispersal_table$pop == .x)) %>%
     purrr::map_int(~ .x[1])
   dispersal_table[indices, ]$tdispersal_gen <- 1

R/interface.R

@@ -13,8 +13,8 @@
 #' @param name Name of the population
 #' @param time Time of the population's first appearance
 #' @param N Number of individuals at the time of first appearance
-#' @param parent Parent population object or "ancestor" (indicating that the
-#'   population does not have an ancestor, and that it is the first population
+#' @param parent Parent population object or \code{NULL} (which indicates that
+#'   the population does not have an ancestor, as it is the first population
 #'   in its "lineage")
 #' @param map Object of the type \code{slendr_map} which defines the world
 #'   context (created using the \code{world} function). If the value
@@ -47,7 +47,7 @@
 #' @export
 #'
 #' @example man/examples/model_definition.R
-population <- function(name, time, N, parent = "ancestor", map = FALSE,
+population <- function(name, time, N, parent = NULL, map = FALSE,
                        center = NULL, radius = NULL, polygon = NULL,
                        remove = NULL, intersect = TRUE,
                        competition = NA, mating = NA, dispersal = NA,
@@ -56,7 +56,7 @@ population <- function(name, time, N, parent = "ancestor", map = FALSE,
     stop("A population name must be a character scalar value", call. = FALSE)
 
   # is this the first population defined in the model?
-  if (is.character(parent) && parent == "ancestor") {
+  if (is.null(parent)) {
     if (!is.logical(map) && !inherits(map, "slendr_map"))
       stop("Ancestral population must specify its 'map'", call. = FALSE)
     else
@@ -95,8 +95,8 @@ population <- function(name, time, N, parent = "ancestor", map = FALSE,
   # keep a record of the parent population
   if (inherits(parent, "slendr_pop"))
     attr(pop, "parent") <- parent
-  else if (is.character(parent) & parent == "ancestor")
-    attr(pop, "parent") <- "ancestor"
+  else if (is.null(parent))
+    attr(pop, "parent") <- "__pop_is_ancestor"
   else
     stop("Suspicious parent population", call. = FALSE)
 

R/print.R

@@ -271,7 +271,7 @@ print_pop_history <- function(x) {
       cat("  - time ")
       cat(event$time, ": ", sep = "")
       parent <- attr(x, "parent")
-      if (is.character(parent) && parent == "ancestor")
+      if (is.character(parent) && parent == "__pop_is_ancestor")
         cat("created as an ancestral population", sprintf("(N = %d)", event$N))
       else {
         cat("split from", parent$pop[1], sprintf("(N = %d)", event$N))

R/utils.R

@@ -601,7 +601,7 @@ check_location_bounds <- function(locations, map) {
 # or ancestry proportions over time) back to user-specified time units
 # (either forward or backward)
 convert_slim_time <- function(times, model) {
-  ancestors <- dplyr::filter(model$splits, parent == "ancestor")
+  ancestors <- dplyr::filter(model$splits, parent == "__pop_is_ancestor")
 
   if (model$direction == "backward") {
     oldest_time <- max(ancestors[, ]$tsplit_orig)

R/visualization.R

@@ -516,7 +516,7 @@ animate_model <- function(model, file, steps, gif = NULL, width = 800, height =
 
 # Create a table of population split edges for graph visualization
 get_split_edges <- function(split_table) {
-  split_edges <- split_table[split_table$parent != "ancestor",
+  split_edges <- split_table[split_table$parent != "__pop_is_ancestor",
                              c("parent", "pop", "tsplit")]
   names(split_edges) <- c("from", "to", "time")
 

README.Rmd

@@ -162,8 +162,7 @@ Note that in order to make this example executable in a reasonable time on my ex
 
 ```{r, message = FALSE}
 afr <- population( # African ancestral population
-  "AFR", parent = "ancestor", time = 52000, N = 3000,
-  map = map, polygon = africa
+  "AFR", time = 52000, N = 3000, map = map, polygon = africa
 )
 
 ooa <- population( # population of the first migrants out of Africa

README.md

@@ -157,7 +157,7 @@ Note that in order to make this example executable in a reasonable time on my ex
 
 ```r
 afr <- population( # African ancestral population
-  "AFR", parent = "ancestor", time = 52000, N = 3000,
+  "AFR", time = 52000, N = 3000,
   map = map, polygon = africa
 )
 

inst/extdata/models/introgression/script.py

@@ -31,7 +31,7 @@ def simulate(
   """Run a slendr simulation and return a tree sequence."""
 
   logging.info("Setting up populations")
-  
+
   # set up demographic history
   demography = msprime.Demography()
   for pop in populations.itertuples():
@@ -44,7 +44,7 @@ def simulate(
           initial_size = resize_events.tail(1).N.values[0]
       else:
           initial_size = pop.N
-  
+
       logging.info(f"Setting up population {pop.pop} with Ne {initial_size}")
       demography.add_population(
           name=pop.pop,
@@ -54,13 +54,13 @@ def simulate(
       # for non-ancestral populations, specify the correct split event and re-set
       # the effective population size (by default in msprime inherited from the
       # parent population)
-      if pop.parent != "ancestor":
+      if pop.parent != "__pop_is_ancestor":
           demography.add_population_split(
               time=length - pop.tsplit_gen,
               derived=[pop.pop],
               ancestral=pop.parent
           )
-  
+
   if not isinstance(samples, pandas.DataFrame) or not len(samples):
       logging.info("No sampling schedule given, generating one automatically")
       samples = pandas.DataFrame(
@@ -70,16 +70,16 @@ def simulate(
           ],
           columns=["n", "pop", "time_gen", "x", "y", "time_orig", "x_orig", "y_orig"]
       )
-  
+
   logging.info("Loading the sampling schedule")
   sample_set = [
       msprime.SampleSet(
          row.n, population=row.pop, time=length - row.time_gen + 1, ploidy=2
       ) for row in samples.itertuples(index=False)
   ]
-  
+
   logging.info("Setting up population resize events")
-  
+
   # schedule population size changes
   for event in resizes.itertuples(index=False):
       if event.how == "step":
@@ -107,9 +107,9 @@ def simulate(
           )
       else:
           sys.exit(f"Unknown event type '{event.how}'")
-  
+
   logging.info("Setting up gene flow events")
-  
+
   # schedule gene flow events
   for event in geneflows.itertuples():
       tstart = length - event.tend_gen + 1
@@ -127,24 +127,24 @@ def simulate(
           source=event.to,
           dest=event._1,
       )
-  
+
   # make sure all slendr events are sorted by time of occurence
   # (otherwise msprime complains)
   demography.sort_events()
-  
+
   if debug:
       print(demography.debug())
-  
+
   logging.info("Running the simulation and generating a tree sequence")
-  
+
   ts = msprime.sim_ancestry(
       samples=sample_set,
       demography=demography,
       sequence_length=sequence_length,
       recombination_rate=recombination_rate,
       random_seed=seed
   )
-  
+
   # compile a set of slendr metadata to be stored in the tree sequence
   slendr_metadata = {
       "slendr": {
@@ -169,15 +169,15 @@ def simulate(
           }
       }
   }
-  
+
   tables = ts.dump_tables()
   tables.metadata_schema = tskit.MetadataSchema({"codec": "json"})
   tables.metadata = slendr_metadata
-  
+
   ts_metadata = tables.tree_sequence()
-  
+
   logging.info("DONE")
-  
+
   return ts_metadata
 
 
@@ -202,48 +202,48 @@ def simulate(
                       help="Print detailed logging information?")
   parser.add_argument("--debug", action="store_true", default=False,
                       help="Print detailed debugging information?")
-  
+
   args = parser.parse_args()
-  
+
   if args.verbose:
       logging.basicConfig(level=logging.INFO)
-  
+
   model_dir = os.path.expanduser(args.model)
-  
+
   logging.info(f"Loading slendr model configuration files from {model_dir}")
-  
+
   if not args.output:
       args.output = pathlib.Path(model_dir, "output_msprime_ts.trees")
-  
+
   if not os.path.exists(model_dir):
       sys.exit(f"Model directory {model_dir} does not exist")
-  
+
   # paths to mandatory slendr configuration files
   populations_path = pathlib.Path(model_dir, "populations.tsv")
   resizes_path = pathlib.Path(model_dir, "resizes.tsv")
   geneflows_path = pathlib.Path(model_dir, "geneflow.tsv")
   length_path = pathlib.Path(model_dir, "length.txt")
   direction_path = pathlib.Path(model_dir, "direction.txt")
   description_path = pathlib.Path(model_dir, "description.txt")
-  
+
   # read model configuration files
   populations = pandas.read_table(populations_path)
   resizes = pandas.DataFrame()
   geneflows = pandas.DataFrame()
   samples = pandas.DataFrame()
-  
+
   if os.path.exists(resizes_path):
       resizes = pandas.read_table(resizes_path)
   if os.path.exists(geneflows_path):
       geneflows = pandas.read_table(geneflows_path)
-  
+
   # read the total simulation length
   length = int(float(open(length_path, "r").readline().rstrip()))
-  
+
   # read the direction of the model ("forward" or "backward")
   direction = open(direction_path, "r").readline().rstrip()
   logging.info(f"Loaded model is specified in the {direction} direction")
-  
+
   # read the description of the model
   description = open(description_path, "r").readline().rstrip()
   logging.info(f"Model description: {description}")
@@ -255,7 +255,7 @@ def simulate(
   ts = simulate(args.sequence_length, args.recombination_rate, args.seed,
                 populations, resizes, geneflows, length, direction, description,
                 samples, args.debug)
-  
+
   output_path = os.path.expanduser(args.output)
 
   logging.info(f"Saving tree sequence output to {output_path}")

inst/extdata/models/introgression/script.slim

@@ -538,7 +538,7 @@ function (void) required_arg(s arg) {
 function (void) create_pop(object<Dictionary> pop) {
     pop_id = pop.getValue("pop_id");
 
-    if (pop.getValue("parent") ==  "ancestor") {
+    if (pop.getValue("parent") ==  "__pop_is_ancestor") {
         log_output("creating " + pop.getValue("pop") + "(p" + pop.getValue("pop_id") + ")");
         sim.addSubpop(pop_id, pop.getValue("N"));
     } else {