Make it easier to run SLiM and msprime on the command line (close #127)

bodkan · Apr 26, 2023 · 2e5b857 · 2e5b857
1 parent e0cf927
commit 2e5b857
Show file tree

Hide file tree

Showing 10 changed files with 552 additions and 455 deletions.
diff --git a/R/cli-runners.R b/R/cli-runners.R
diff --git a/R/compilation.R b/R/compilation.R
diff --git a/R/engines.R b/R/engines.R
diff --git a/docs/reference/msprime.html b/docs/reference/msprime.html
diff --git a/docs/reference/slim.html b/docs/reference/slim.html
diff --git a/man/msprime.Rd b/man/msprime.Rd
diff --git a/man/slim.Rd b/man/slim.Rd
diff --git a/tests/testthat/test-engines.R b/tests/testthat/test-engines.R
@@ -0,0 +1,92 @@
+test_that("only serialized models can be run on the command line", {
+  pop1 <- population("pop1", N = 1000, time = 1)
+  pop2 <- population("pop2", N = 1000, time = 2, parent = pop1)
+
+  model <- compile_model(list(pop1, pop2), generation_time = 1, direction = "forward", simulation_length = 1000, serialize = FALSE)
+
+  out <- tempfile()
+  expect_error(
+    msprime(model, sequence_length = 1e6, recombination_rate = 1e-8, output = out, run = FALSE),
+    "Impossible to run a non-serialized slendr model on the command line"
+  )
+})
+
+test_that("msprime command run manually on the command line give the correct output", {
+  pop1 <- population("pop1", N = 1000, time = 1)
+  pop2 <- population("pop2", N = 1000, time = 2, parent = pop1)
+
+  model <- compile_model(list(pop1, pop2), generation_time = 1, direction = "forward", simulation_length = 1000)
+
+  # check that a simulated tree-sequence file is where it's supposed to be if the
+  # model is run on the CLI
+  out_cmd <- capture.output(
+    cmd <- msprime(model, sequence_length = 1e6, recombination_rate = 1e-8, output = out, run = FALSE, random_seed = 42))
+  system(cmd, ignore.stdout = TRUE)
+  ts_manual <- ts_load(out, model = model)
+  expect_s3_class(ts_manual, "slendr_ts")
+
+  # check that the manually simulated tree-sequence matches what comes from running inside slendr
+  ts_r <- msprime(model, sequence_length = 1e6, recombination_rate = 1e-8, random_seed = 42)
+  expect_equal(ts_nodes(ts_manual), ts_nodes(ts_r))
+
+  # and check that the command itself matches
+  expect_equal(gsub(" *$", "", out_cmd[4]), gsub(" *$", "", cmd))
+})
+
+test_that("slim command run manually on the command line give the correct output", {
+  pop1 <- population("pop1", N = 1000, time = 1)
+  pop2 <- population("pop2", N = 1000, time = 2, parent = pop1)
+
+  model <- compile_model(list(pop1, pop2), generation_time = 1, direction = "forward", simulation_length = 1000)
+
+  # check that a simulated tree-sequence file is where it's supposed to be if the
+  # model is run on the CLI
+  out_cmd <- capture.output(
+    cmd <- slim(model, sequence_length = 1e6, recombination_rate = 1e-8, output = out, run = FALSE, random_seed = 42))
+  system(cmd, ignore.stdout = TRUE)
+  ts_manual <- ts_load(out, model = model)
+  expect_s3_class(ts_manual, "slendr_ts")
+
+  # check that the manually simulated tree-sequence matches what comes from running inside slendr
+  ts_r <- slim(model, sequence_length = 1e6, recombination_rate = 1e-8, random_seed = 42)
+  expect_equal(ts_nodes(ts_manual), ts_nodes(ts_r))
+
+  # and check that the command itself matches
+  expect_equal(
+    gsub(" *$", "", out_cmd[-c(1:3, length(out_cmd) - 1, length(out_cmd))]),
+    strsplit(cmd, "\n")[[1]]
+  )
+})
+
+test_that("ensure that a model reaches full coalescence", {
+  pop1 <- population("pop1", N = 1000, time = 1)
+  pop2 <- population("pop2", N = 1000, time = 1)
+  model <- compile_model(list(pop1, pop2), generation_time = 1, direction = "forward", simulation_length = 1000)
+
+  expect_error(msprime(model, sequence_length = 1e6, recombination_rate = 1e-8, verbose = TRUE),
+               "Multiple ancestral populations without a common ancestor")
+
+  pop3 <- population("pop2", N = 1000, time = 2, parent = pop1)
+  model <- compile_model(list(pop1, pop3), generation_time = 1, direction = "forward", simulation_length = 1000)
+
+  expect_s3_class(
+    msprime(model, sequence_length = 1e6, recombination_rate = 1e-8, verbose = FALSE),
+    "slendr_ts"
+  )
+})
+
+test_that("population size must be a non-negative number", {
+  expect_error(pop <- population("asd", N = -1, time = 100),
+               "Population size must be a non-negative number")
+  expect_error(pop <- population("asd", N = 0, time = 100),
+               "Population size must be a non-negative number")
+  expect_s3_class(pop <- population("asd", N = 1, time = 100), "slendr_pop")
+})
+
+test_that("population time must be a non-negative number", {
+  expect_error(pop <- population("asd", N = 100, time = -1),
+               "Split time must be a non-negative number")
+  expect_error(pop <- population("asd", N = 100, time = 0),
+               "Split time must be a non-negative number")
+  expect_s3_class(pop <- population("asd", N = 1, time = 100), "slendr_pop")
+})
diff --git a/tests/testthat/test-msprime.R b/tests/testthat/test-msprime.R
@@ -323,36 +323,3 @@ test_that("AFS distributions from SLiM and msprime simulations match", {
 
   expect_equal(afs, orig_afs)
 })
-
-test_that("ensure that a model reaches full coalescence", {
-  pop1 <- population("pop1", N = 1000, time = 1)
-  pop2 <- population("pop2", N = 1000, time = 1)
-  model <- compile_model(list(pop1, pop2), generation_time = 1, direction = "forward", simulation_length = 1000)
-
-  expect_error(msprime(model, sequence_length = 1e6, recombination_rate = 1e-8, verbose = TRUE),
-               "Multiple ancestral populations without a common ancestor")
-
-  pop3 <- population("pop2", N = 1000, time = 2, parent = pop1)
-  model <- compile_model(list(pop1, pop3), generation_time = 1, direction = "forward", simulation_length = 1000)
-
-  expect_s3_class(
-    msprime(model, sequence_length = 1e6, recombination_rate = 1e-8, verbose = FALSE),
-    "slendr_ts"
-  )
-})
-
-test_that("population size must be a non-negative number", {
-  expect_error(pop <- population("asd", N = -1, time = 100),
-               "Population size must be a non-negative number")
-  expect_error(pop <- population("asd", N = 0, time = 100),
-               "Population size must be a non-negative number")
-  expect_s3_class(pop <- population("asd", N = 1, time = 100), "slendr_pop")
-})
-
-test_that("population time must be a non-negative number", {
-  expect_error(pop <- population("asd", N = 100, time = -1),
-               "Split time must be a non-negative number")
-  expect_error(pop <- population("asd", N = 100, time = 0),
-               "Split time must be a non-negative number")
-  expect_s3_class(pop <- population("asd", N = 1, time = 100), "slendr_pop")
-})
diff --git a/vignettes/vignette-04-nonspatial-models.Rmd b/vignettes/vignette-04-nonspatial-models.Rmd
@@ -65,8 +65,7 @@ The compilation step is also the same. The only (internal) difference is that we
 ```{r}
 model <- compile_model(
   populations = list(afr, ooa, ehg, eur, ana, yam),
-  gene_flow = gf, generation_time = 30,
-  path = paste0(tempfile(), "_non-spatial-model")
+  gene_flow = gf, generation_time = 30
 )
 ```