Add option to name time units used by the model

R/compilation.R

@@ -30,6 +30,7 @@
 #'   present? Useful for non-interactive uses. In an interactive mode, the user
 #'   is asked to confirm the deletion manually.
 #' @param description Optional short description of the model
+#' @param time_units Units of time in which model event times are to be interpreted
 #' @param resolution How many distance units per pixel?
 #' @param competition,mating Maximum spatial competition and mating choice
 #'   distance
@@ -54,7 +55,7 @@ compile_model <- function(
     populations, generation_time, gene_flow = list(),
     direction = NULL, simulation_length = NULL,
     serialize = TRUE, path = NULL, overwrite = FALSE, force = FALSE,
-    description = "",
+    description = "", time_units = "",
     resolution = NULL, competition = NULL, mating = NULL, dispersal = NULL,
     extension = NULL
 ) {
@@ -232,7 +233,7 @@ setting `direction = 'backward'.`", call. = FALSE)
     checksums <- write_model(
       path, populations, admix_table, map_table, split_table, resize_table,
       dispersal_table, generation_time, resolution, simulation_length, time_dir, slim_script,
-      description, map
+      description, time_units, map
     )
   else
     checksums <- NULL
@@ -257,6 +258,7 @@ setting `direction = 'backward'.`", call. = FALSE)
     orig_length = simulation_length,
     direction = time_dir,
     description = description,
+    time_units = time_units,
     checksums = checksums,
     customized = customized
   )
@@ -303,6 +305,7 @@ read_model <- function(path) {
   path_orig_length <- file.path(path, "orig_length.txt")
   path_direction <- file.path(path, "direction.txt")
   path_description <- file.path(path, "description.txt")
+  path_time_units <- file.path(path, "time_units.txt")
 
   if (!dir.exists(path))
     stop(sprintf("Model directory '%s' does not exist", path), call. = FALSE)
@@ -314,7 +317,8 @@ read_model <- function(path) {
   generation_time <- scan(path_generation_time, what = integer(), quiet = TRUE)
   length <- as.integer(scan(path_length, what = numeric(), quiet = TRUE))
   orig_length <- as.integer(scan(path_orig_length, what = numeric(), quiet = TRUE))
-  description <- scan(path_description, what = character(), quiet = TRUE)
+  description <- readLines(path_description)
+  time_units <- readLines(path_time_units)
 
   split_table <- utils::read.table(path_splits, header = TRUE, stringsAsFactors = FALSE)
   resize_table <- NULL
@@ -357,6 +361,8 @@ read_model <- function(path) {
     length = length,
     orig_length = orig_length,
     direction = direction,
+    description = description,
+    time_units = time_units,
     checksums = checksums,
     customized = customized
   )
@@ -399,7 +405,7 @@ verify_checksums <- function(files, hashes) {
 write_model <- function(path, populations, admix_table, map_table, split_table,
                         resize_table, dispersal_table,
                         generation_time, resolution, length, direction,
-                        script_source, description, map) {
+                        script_source, description, time_units, map) {
   saved_files <- c()
 
   # table of split times and initial population sizes
@@ -457,11 +463,13 @@ write_model <- function(path, populations, admix_table, map_table, split_table,
   saved_files["orig_length"] <- file.path(path, "orig_length.txt")
   saved_files["direction"] <- file.path(path, "direction.txt")
   saved_files["description"] <- file.path(path, "description.txt")
+  saved_files["time_units"] <- file.path(path, "time_units.txt")
   base::write(generation_time, file.path(path, "generation_time.txt"))
   base::write(round(length / generation_time), file.path(path, "length.txt"))
   base::write(length, file.path(path, "orig_length.txt"))
   base::write(direction, file.path(path, "direction.txt"))
   base::write(description, file.path(path, "description.txt"))
+  base::write(time_units, file.path(path, "time_units.txt"))
 
   saved_files["slim_script"] <- file.path(path, "script.slim")
   saved_files["msprime_script"] <- file.path(path, "script.py")

R/visualization.R

@@ -448,10 +448,8 @@ plot_model <- function(model, sizes = TRUE, proportions = FALSE, gene_flow = TRU
   } else
     gene_flow <- NULL
 
-  if (model$direction == "forward")
-    ylabel <- "time since the start"
-  else
-    ylabel <- "time before present"
+  ylabel <- "time"
+  if (model$time_units != "") ylabel <- sprintf("%s (%s)", ylabel, model$time_units)
 
   if (log) ylabel <- paste(ylabel, "(log scale)")
 

tests/testthat/test-serialization.R

@@ -7,12 +7,13 @@ test_that("read_model() restores a single-map model object", {
 
   model_dir <- file.path(tempdir(), "tmp-single-map-model-serialization")
   model1 <- compile_model(pop, path = model_dir, resolution = 10000, generation_time = 1, overwrite = TRUE, force = TRUE,
-                    competition = 100e3, mating = 100e3, dispersal = 10e3, direction = "backward")
+                    competition = 100e3, mating = 100e3, dispersal = 10e3, direction = "backward",
+                    description = "This is a model description", time_units = "These are the model time units")
   model2 <- read_model(model1$path)
 
   # make sure that all components of the model list object before and after
   # serialization are equal
-  components <- c("checksums", "splits", "resizes", "geneflow", "maps", "direction", "length", "generation_time", "resolution", "world")
+  components <- c("checksums", "splits", "resizes", "geneflow", "maps", "direction", "length", "generation_time", "resolution", "world", "description", "time_units")
   expect_true(all(sapply(components, function(i) all.equal(model1[[i]], model2[[i]]))))
   expect_true(all(sapply(seq_along(model1$populations), function(i) all(model1$populations[[i]] == model2$populations[[i]]))))
   expect_equal(names(model1$populations), names(model1$populations))
@@ -47,7 +48,7 @@ test_that("read_model() restores a complex model object", {
   )
   model2 <- read_model(model1$path)
 
-  components <- c("checksums", "splits", "resizes", "geneflow", "maps", "length", "orig_length", "direction", "generation_time", "resolution", "world")
+  components <- c("checksums", "splits", "resizes", "geneflow", "maps", "length", "orig_length", "direction", "generation_time", "resolution", "world", "description", "time_units")
   expect_true(all(sapply(components, function(i) all.equal(model1[[i]], model2[[i]]))))
   expect_true(all(sapply(seq_along(model1$populations), function(i) all(model1$populations[[i]] == model2$populations[[i]]))))
   expect_equal(names(model1$populations), names(model1$populations))
@@ -83,7 +84,7 @@ test_that("read_model() restores a single-map model object (nonspatial)", {
 
   # make sure that all components of the model list object before and after
   # serialization are equal
-  components <- c("checksums", "splits", "resizes", "geneflow", "maps", "direction", "length", "generation_time", "resolution", "world")
+  components <- c("checksums", "splits", "resizes", "geneflow", "maps", "direction", "length", "generation_time", "resolution", "world", "description", "time_units")
   expect_true(all(sapply(components, function(i) all.equal(model1[[i]], model2[[i]]))))
   expect_true(all(sapply(seq_along(model1$populations), function(i) all(unlist(model1$populations[[i]]) == unlist(model2$populations[[i]])))))
 })
@@ -106,11 +107,13 @@ test_that("read_model() restores a complex model object (nonspatial)", {
     populations = list(p1, p2, p3, p4, p5),
     gene_flow = geneflow,
     generation_time = 30,
-    overwrite = TRUE, force = TRUE
+    overwrite = TRUE, force = TRUE,
+    description = "This is a description",
+    time_units = "These are time units"
   )
   model2 <- read_model(model1$path)
 
-  components <- c("checksums", "splits", "resizes", "geneflow", "maps", "length", "orig_length", "direction", "generation_time", "resolution", "world")
+  components <- c("checksums", "splits", "resizes", "geneflow", "maps", "length", "orig_length", "direction", "generation_time", "resolution", "world", "description", "time_units")
   expect_true(all(sapply(components, function(i) all.equal(model1[[i]], model2[[i]]))))
   expect_true(all(sapply(seq_along(model1$populations), function(i) all(unlist(model1$populations[[i]]) == unlist(model2$populations[[i]])))))
   expect_equal(names(model1$populations), names(model1$populations))