search

bodkan / slendr

Population genetic simulations in R 🌍
https://bodkan.net/slendr
Other
55 stars 6 forks source link

Locking population resizes with expansions/contractions interferes with determining the direction of time #143

Closed bodkan closed 11 months ago

bodkan commented 1 year ago

Consider this chunk of code which determines the spatio-temporal dynamics of population expansion:

ana <- population( # Anatolian farmers
  name = "ANA", time = 13000, N = 5000, map = map,
  center = c(34, 38), radius = 500e3, polygon = anatolia
) %>%
  expand_range(
    by = 3000e3, start = 10000, end = 7000,
    polygon = join(europe, anatolia), snapshots = 20
  )

This creates the following object -- no issues there:

> ana
slendr 'population' object 
-------------------------- 
name: ANA 
habitat: terrestrial

number of spatial maps: 21 
map: internal coordinate reference system EPSG 3035 
stays until the end of the simulation

population history overview:
  - time 13000: created as an ancestral population (N = 5000)
  - time 10000-7000: range expansion

However, when lock = TRUE is set in expand_range() (an option which allows automatic adjustment of population sizes to maintain +/- same densities), this causes an issue with the direction of times of individual events:

ana <- population( # Anatolian farmers
  name = "ANA", time = 13000, N = 5000, map = map,
  center = c(34, 38), radius = 500e3, polygon = anatolia
) %>%
  expand_range(
    by = 3000e3, start = 10000, end = 7000,
    polygon = join(europe, anatolia), snapshots = 20,
    lock = TRUE  # <---- this right here
  )
> ana
slendr 'population' object 
-------------------------- 
name: ANA 
habitat: terrestrial

number of spatial maps: 21 
map: internal coordinate reference system EPSG 3035 
stays until the end of the simulation

population history overview:
  - time 13000: created as an ancestral population (N = 5000)
  - time 10000-7000: range expansion
  - time 9850: resize from 5000 to 5960 individuals
  - time 9700: resize from 5960 to 7278 individuals
  - time 9550: resize from 7278 to 8743 individuals
  - time 9400: resize from 8743 to 10282 individuals
  - time 9250: resize from 10282 to 11849 individuals
  - time 9100: resize from 11849 to 13451 individuals
  - time 8950: resize from 13451 to 15103 individuals
  - time 8800: resize from 15103 to 16810 individuals
  - time 8650: resize from 16810 to 18574 individuals
  - time 8500: resize from 18574 to 20425 individuals
  - time 8350: resize from 20425 to 22380 individuals
  - time 8200: resize from 22380 to 24438 individuals
  - time 8050: resize from 24438 to 26558 individuals
  - time 7900: resize from 26558 to 28690 individuals
  - time 7750: resize from 28690 to 30839 individuals
  - time 7600: resize from 30839 to 33008 individuals
  - time 7450: resize from 33008 to 35199 individuals
  - time 7300: resize from 35199 to 37414 individuals
  - time 7150: resize from 37414 to 39662 individuals
  - time 7000: resize from 39662 to 41898 individuals

Setting lock = TRUE simply auto-schedules normal resize() events, which seems fine, except that the orders of times of events is then no longer monotonic: 13000, 10000, 7000, 9850, 9700...

This shouldn't be a problem by itself, except that slendr's internal checks of time direction require that times of subsequent events are monotonic. That's how we can let the users specify time units in whatever direction they want -- slendr can always translate to forward-ticks SLiM time or backwards-generations msprime time in a "smart way".

This should be fixable, but it will probably make the time-direction checking code even more complex. 🙄

This honestly isn't a backbreaking bug, just a small annoyance.

For completeness, here are the spatial map components needed for the code above to work:

map <- world(xrange = c(-13, 70), yrange = c(18, 65), crs = "EPSG:3035")

europe <- region(
  "Europe", map,
  polygon = list(
    c(-8, 35), c(-5, 36), c(10, 38), c(20, 35), c(25, 35),
    c(33, 45), c(20, 58), c(-5, 60), c(-15, 50)
  )
)

anatolia <- region(
  "Anatolia", map,
  polygon = list(c(28, 35), c(40, 35), c(42, 40),
                 c(30, 43), c(27, 40), c(25, 38))
)
bodkan commented 11 months ago

So I think I fixed it?

In the end it doesn't seem like dramatic changes to the internals were needed. Turns out that if an range expansion/shrinking is scheduled with locked size change, the size change progression can be used to determine the model time direction, skipping the apparent inconsistency caused by the range change and population size events. What I did was to simply skip the range expansion event from the loop used to determine time direction if the range change is associated with population size change. If it isn't, then the code works just as before. The change in the code can be found here.

The tests and examples run as before, so hopefully I didn't break anything else in the process.