---
title: "Traditional, non-spatial models"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteEncoding{UTF-8}
  %\VignetteIndexEntry{Traditional, non-spatial models}
  %\VignetteEngine{knitr::rmarkdown}
---

```{r, include = FALSE}
env_present <- slendr:::is_slendr_env_present()

knitr::opts_chunk$set(
  collapse = FALSE,
  comment = "#>",
  eval = env_present,
  fig.width = 6,
  fig.height = 4,
  dpi = 70
)
```

The biggest selling point of the *slendr* package is that you can program spatiotemporal population genetics models in R and have them execute automatically in SLiM. However, there are several reasons why you might be interested in using *slendr* even for non-spatial models. First, R is a language that many scientists already know, and being able to simulate data from the comfort of an R interface significantly lowers the barrier of entry. Second, because *slendr* makes SLiM appear almost as if it were just another R library, running simulations (spatial *and* non-spatial), fitting models, exploring parameter grids, calculating statistics, and visualization of results can be all be performed without leaving the R interface.

In this vignette, we will demonstrate how to program non-spatial models in *slendr*. However, we should start by noting that there is almost no difference between code for non-spatial and spatial models in *slendr*. The only visible difference is that spatial models include a `map =` argument in the `population()` constructor function of ancestral population(s), and non-spatial models do not. That's it, that's the difference. Switching between spatial and non-spatial models is performed internally by the package, without any user intervention.

To make the comparison clearer, we will use the example from the *slendr* [landing page](https://github.com/bodkan/slendr), but we will implement it in a non-spatial context (i.e., as a traditional random mating simulation).

First, let's define population objects, splits, and other demographic events (note the missing `map` argument, which is set to `FALSE` by default):

```{r}
library(slendr)

init_env()

# African ancestral population
afr <- population("AFR", time = 100000, N = 3000)

# first migrants out of Africa
ooa <- population("OOA", parent = afr, time = 60000, N = 500, remove = 23000) %>%
  resize(N = 2000, time = 40000, how = "step")

# Eastern hunter-gatherers
ehg <- population("EHG", parent = ooa, time = 28000, N = 1000, remove = 6000)

# European population
eur <- population("EUR", parent = ehg, time = 25000, N = 2000) %>%
  resize(N = 10000, how = "exponential", time = 5000, end = 0)

# Anatolian farmers
ana <- population("ANA", time = 28000, N = 3000, parent = ooa, remove = 4000)

# Yamnaya steppe population
yam <- population("YAM", time = 7000, N = 500, parent = ehg, remove = 2500)
```

We can define gene flow events in the same way as we did for the spatial model:

```{r}
gf <- list(
  gene_flow(from = ana, to = yam, rate = 0.5, start = 6500, end = 6400),
  gene_flow(from = ana, to = eur, rate = 0.5, start = 8000, end = 6000),
  gene_flow(from = yam, to = eur, rate = 0.75, start = 4000, end = 3000)
)
```

The compilation step is also the same. The only (internal) difference is that we skip the rasterization of vector maps that is performed for spatial models in order to control and restrict population boundaries:

```{r}
model <- compile_model(
  populations = list(afr, ooa, ehg, eur, ana, yam),
  gene_flow = gf, generation_time = 30,
  time_units = "years before present"
)
```

Using the `plot_map()` function doesn't make sense, as there are no spatial maps to plot. However, we can still plot the demographic graph, verifying that the model has been specified correctly using the function `plot_model()` as shown below.

Let's say we also want to schedule specific sampling events at specific times (which record only specified individuals in a tree sequence). We can use `schedule_sampling()` to do just that:

```{r}
samples <- schedule_sampling(
  model,
  times = c(0, 5000, 12000, 20000, 35000, 39000, 43000),
  list(eur, 3), list(ehg, 1), list(yam, 1), list(ana, 3), list(ooa, 1), list(afr, 1)
)
```

Because it's not just the model itself that's useful to visually verify (which is the main purpose of `plot_model()`) but also the sampling scheme, _slendr_ makes it possible to overlay the numbers of individuals scheduled for tree-sequence recording from each lineage at each timepoint:

```{r, non-spatial-graph_sampling, fig.width=6, fig.height=4}
plot_model(model, samples = samples)
```

Even the final step---execution of the model in SLiM---is the same, using the built-in `slim()` function:

```{r, eval = FALSE}
ts_slim <- slim(model, sequence_length = 100000, recombination_rate = 0)
```

Even for non-spatial models, this function still uses the same SLiM back end script used for spatial models. The only difference is that all spatial features are switched off, making the model run as a simple random-mating simulation.

Given that we are running a non-spatial simulation, you might wonder if it wouldn't be more efficient to use a coalescent simulator. Indeed, *slendr* also provides an alternative *msprime* back end just for this purpose. We could run the exact same simulation with *msprime* like this:

```{r}
ts_msprime <- msprime(model, sequence_length = 100000, recombination_rate = 0)
```

In fact, because both SLiM and *msprime* back ends save outputs in a tree sequence format, we can analyse them using the same tools. See [this vignette](https://www.slendr.net/articles/vignette-05-tree-sequences.html) for more information about tree sequence analysis with *slendr*, and for more discussion on alternative simulation back ends and more extensive examples of data analysis with tree sequences you can read [this tutorial](https://www.slendr.net/articles/vignette-07-backends.html).

## Extracting parameters from a model or tree sequences

In some situations (such as when model parameters are drawn from random distributions and we need to know which parameters had been used _after_ the simulations), function `extract_parameters()` can be used. This function peeks into a _slendr_ tree sequence object and extract parameters of the original _slendr_ model:

```{r}
extract_parameters(ts_msprime)
```

For completeness, although this isn't relevant in this example either, the function can also get parameters from any compiled model. For instance, we can check which parameters were used to compile the built-in _slendr_ introgression model by running:

```{r}
introgression_model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

extract_parameters(introgression_model)
```

As we can see, `extract_parameters()` returns a list of data frames, one data frame for each aspect of a demographic model (where applicable).