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}
end
```

We can check that we are able to create a `Torus{Bool}`

with 300 × 150 cells:

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

```
Main.WeaveSandBox2.Torus{Bool}
```

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 162 methods)
```

Note the definition of `eltype`

- Julia’s type system is really cool. We can
next check that `size`

returns the correct value:

```
size(world)
```

```
(150, 300)
```

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

```
world_update = similar(world)
size(world_update)
```

```
(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...]
end
```

```
getindex (generic function with 378 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)
end
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...)
end
```

```
setindex! (generic function with 154 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]
end
```

```
getindex (generic function with 379 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
end
occupancy(world)
```

```
0.09913333333333334
```

At this point, our updated world is still empty:

```
occupancy(world_update)
```

```
0.0
```

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
else
# Cell empty - becames empty
world_update[i,j] = rand() < P
end
end
end
# When the generation is done, we update the world
for i in eachindex(world)
world[i] = world_update[i]
world_update[i] = false
end
end
```

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:

```
occupancy(world)
```

```
0.12935555555555556
```

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.