Let’s play a game. We will write a model of population dynamics, and repeat it over a number of timesteps; then we will change the parameters, and draw a bifurcation diagram. These is the goal of the game. The only rule is as follows: if we write a loop, we lose. That’s right, we are going to do population dynamics using only functional programming. Because it is fun.
Setting up the model(s)
Let’s begin by defining two functions: one for exponential growth, and one for regulation. These are function generators, in that they will return a function, which can then be called with parameters.
R(r::Float64) = (n) -> n*r K(k::Float64) = (n) -> n*(1.0-n/k)
We can create a growth model from these two components, by calling
R. This is strictly useless, but a good opportunity to showcase the
\circ) operator, which in Julia works exactly like it works in the
mathematical world, i.e. it composes functions.
growth = (R(2.3)∘K(1.0))
Let’s check what this does. First, it calls
K to generate a function,
which is $g(n) = 1-n\times K$, where $K$ is the carrying capacity, and $n$
is the population density. Then, it calls
R to generate a function $f(n) =
r\times Y$, where $Y$ is actually $g(n)$. So the final function is $f\circ g(n)
= f(g(n))$, which is $r\times n\times(1 - n/K)$. This is the logistic model.
We can check the result of
growth(1.0) (0.0), and
(0.0), as well as something above the carrying capacity, like
This seems to be working! The next step is to actually run this a number of
times, and to increment the population size. This can be achieved with the
accumulate function. To
accumulate something, we need an operator with
x[i+1]. In our case, we want to add
n, so we will not need the second argument, it’s simply here as a
placeholder. Let us build the correct operator:
model = (current, _) -> current + growth(current)
Let’s apply this to a series of 50 elements, starting with 0.01. As the value of $r$ is large (2.3), this should lead to oscillations with a period of 2:
nt = accumulate(model, zeros(Float64, 50); init=0.01) plot(nt, leg=false, c=:teal, frame=:origin)
Sweet - it works.
Now let’s ramp this up, and see if we can come up with a bifurcaton diagram
without writing a
for loop. To do this, we will need a sequence of values
of $r$ to try, which is always the first step in a bifurcation diagram. As
we have read our textbooks (or because we have used a very elegant approach
to figure it out), we know that the system crashes above $r = 3.0$, so we
will only go up until this point.
r_seq = 2.0:0.001:3.0
Now, we can
map through these values, to create a series of models:
models = map(r -> (current, _) -> current + (R(r)∘K(1.0))(current), r_seq);
Every element in this array is a model, which can then be called within
accumulate, as we did before. At this point, we can generate the bifurcation diagram:
# Generate the timeseries raw_timeseries = map(m -> accumulate(m, zeros(Float64, 5000); init=0.01), models) # Get the unique values across the latest 200 steps endpoints = map(ts -> unique(ts[(end-200:end)]), raw_timeseries) # Generate a list of the values of r for every simulation rs = map(i -> fill(collect(r_seq)[i], length(endpoints[i])), eachindex(endpoints)) # Plot! scatter(rs, endpoints, leg=false, c=:teal, msc=:transparent, msw=0.0, ms=1, alpha=0.2, grid=false, frame=:none)
This is exactly what we were expecting to see!
This is, almost without a doubt, the most needlessly complicated way to solve
this problem. But it goes to show that the tools of functional programming
can be applied even for things that we would traditionally solve with a good
for loop. I particularly like the
\circ notation, which is providing
code that is ever closer to its mathematical meaning.