# Population dynamics and functional programming

## We don't need no iteration

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 `K`

first,
then `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 `growth(0.0)`

(0.0), as well as something above the carrying capacity, like
`growth(1.2)`

(-0.55).

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
two arguments, `x[i]`

and `x[i+1]`

. In our case, we want to add `growth(n)`

to `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.

## Bifurcation diagram

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!

## A conclusion?

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
old `for`

loop. I particularly like the `\circ`

notation, which is providing
code that is ever closer to its mathematical meaning.