Meta-populations on a toroidal world

Using custom types to simplify the code

Last week, we had a class live implementation of a model of metapopulation dynamics, during which we discussed the various ways in which we can deal with the edges of the world. Because time was scarce, we decided to simulate more than we needed, and to focus on the center of the world. In this entry, I would like to highlight ways in which interfaces in Julia can lead to idiomatic code.

Throughout this example, we will work on toroidal worlds. A torus is a grid whose borders are aligned, so that we can wrap around in both directions. Therefore, our typed torus will be a wrapper around a similarly typed array.

mutable struct Torus{T}
  cells::Array{T, 2}

We can check that we are able to create a Torus{Bool} with 300 × 150 cells:

world = Torus(zeros(Bool, (150, 300)))

At this point, we cannot do much with this new object, because we do not have methods that can be applied to a Torus. Therefore, the first step is to define an interface.

Defining the interface

The Julia manual has a very interesting page on Interfaces, which are ways to ensure that user-defined types behave in the same way as the built-in ones. This allows for a lot of composability, and ensures that all well-written code can be used in a Julian way. For our example, we will define the interface for indexing and arrays, and try to avoid the interface for iteration (it would not be more difficult to add it, but we will not need it).

The interface for arrays is very straightforward, as we will work on the cells field as we would on any matrix. We will first define the size, similar, and eltype methods:

import Base.size
size(t::Torus{T}) where {T} = size(t.cells)

import Base.similar
similar(t::Torus{T}) where {T} = Torus(similar(t.cells))

import Base.eltype
eltype(t::Torus{T}) where {T} = T
eltype (generic function with 203 methods)

Note the definition of eltype - Julia’s type system is really cool. We can next check that size returns the correct value:

(150, 300)

At this point, we can also create a similarly sized world for the updates during the simulation:

world_update = similar(world)
(150, 300)

So far, so good. We can now write the indexing interface. The only difference between a torus and a grid is that if the coordinate is lower than 1, we wrap back to the top, and if it is larger than the maximum, we wrap to the bottom. Simply put, this amounts to not using the $x$ position directly, but rather to use $\text{mod}(x-1, M)+1$, where $M$ is the size alongside a dimension. It would be easier to express this in a language where indexing starts at 0, but these don’t make any sense, so I’m happy with needing the +1 and -1.

If we have 200 cells, we can check that this approach works:

X = [0, 201, 200, 1, -201]
M = 200
collect(zip(X, [mod(x-1, M)+1 for x in X]))
5-element Array{Tuple{Int64,Int64},1}:
 (0, 200)
 (201, 1)
 (200, 200)
 (1, 1)
 (-201, 199)

Therefore, we have enough to write the indexing interface. We will simplify the task by assuming that accessing a position using a single integer is not subject to overflow (and we will only use this for cases where it is trivially true anyways):

import Base.getindex
getindex(t::Torus{T}, i::Int) where {T} = t.cells[i]

import Base.eachindex
eachindex(t::Torus{T}) where {T} = eachindex(t.cells)
eachindex (generic function with 27 methods)

We can now write the method using the full position:

function getindex(t::Torus{T}, I::Vararg{Int, N}) where {T, N}
  i = [mod(I[pos]-1, size(t)[pos])+1 for pos in 1:N]
  return t.cells[i...]
getindex (generic function with 487 methods)

The setindex! methods are equally straightforward to write:

import Base.setindex!

function setindex!(t::Torus{T}, v::T, i::Int) where {T}
  setindex!(t.cells, v, i)

function setindex!(t::Torus{T}, v::T, I::Vararg{Int, N}) where {T, N}
  i = [mod(I[pos]-1, size(t)[pos])+1 for pos in 1:N]
  setindex!(t.cells, v, i...)
setindex! (generic function with 232 methods)

At some point, we will also want to select a range of values in our toroidal world (to get the neighbors of a point), so building a method with UnitRange is a good idea:

function getindex(t::Torus{T}, lon::UnitRange{Int}, lat::UnitRange{Int}) where {T}
  ilon = [mod(l-1, size(t)[1])+1 for l in lon]
  ilat = [mod(l-1, size(t)[2])+1 for l in lat]
  return t.cells[ilon, ilat]
getindex (generic function with 488 methods)

Alright! Almost done. Let’s add a function to measure occupancy:

occupancy(t::Torus{T}) where {T} = sum(t.cells .> zero(T))/prod(size(t))
occupancy (generic function with 1 method)

By the way, our function would return the occupancy (proportion of occupied sites) even if our toroidal world had stored numerical values, so this function would work in situations where we model population dynamics, for example.

And now, we are all set!

Starting the simulation

At every timestep, we will look through our world, cell by cell, decide if there is an extinction, and decide if there is a colonisation. We will assume that every occupied neighbor (of which there are $n$) is able to colonize the focal patch at rate $c$, and that the rate of extinction is $\epsilon$. Therefore the probability of a patch being colonized is $P = 1-(1-c)^n$ (the probability of having at least one colonization event is one minus the propability of all colonization attempts failing), from which we get that if a patch is initially occupied, it remains so with probability $(1-\epsilon)+\epsilon P$, and a patch that is initially empty becomes occupied at rate $P$ (this can be checked using a probability tree for the success of colonization and extinction, it’s all good fun). This model assumes that we have colonisation and extinction at the same time, by the way.

We need to decide on an initial occupancy for our world, so let’s say this will be 0.1:

for i in eachindex(world)
  world[i] = rand() < 0.1
  world_update[i] = false

At this point, our updated world is still empty:


We can now perform the actual simulation.

Running the simulation

Once we have picked $c$ and $\epsilon$, the actual simulation loop is trivial. The trick is to work on the probability that a cell becomes occupied as a function of its initial state, which we have calculated before.

c, ε = 0.05, 0.3
for t in 1:100
  for i in 1:size(world)[1]
    for j in 1:size(world)[2]
      # How many neighbors (excluding oneself)
      n = sum(world[(i-1):(i+1),(j-1):(j+1)])-world[i,j]
      # Actual colonization probability
      P = 1-(1-c)^n
      # Choices based on the probability tree
      if world[i,j]
        # Cell occupied - remains occupied
        world_update[i,j] = rand() < (1-ε)+ε*P
        # Cell empty - becames empty
        world_update[i,j] = rand() < P
  # When the generation is done, we update the world
  for i in eachindex(world)
    world[i] = world_update[i]
    world_update[i] = false

And we are done! This particular implementation is very slow (for many reasons), but it works as a rough first draft. We can check the final occupancy:


The only time where we needed the Torus type was here:

n = sum(world[(i-1):(i+1),(j-1):(j+1)])-world[i,j]

When we are near the borders of the world, it is possible to get values of i-1 that are lower than 1, or values of j+1 larger than the maximum. With the Torus type, we are sure that we will get the correct positions, by wrapping around the world. This was, admittedly, overkill, but this was nevertheless a nice illustration of how we can use custom types like base types at the cost of writing a few specialized methods.