N = 1000  # Number of species.
t_end = 100000  # End time of simulation.
t_space = 1:t_end  # An iterable for the time steps. Note that Julia is 1-indexed.

# Create the species' fitness factor.
# rand(N) returns a vector of size N with random elements between 0.0 and 1.0.
species = rand(N)

# Fix boundary conditions
# ir and il are right indices and left indices respectively
# ir[i] gives the right index at position i, ususally i+1.
#
# `collect` converts a range, for example 2:N, into an array.
# It serves the same purpose as doing `list(range(2, N))` in Python.
ir = collect(2:N+1)
il = collect(0:N-1)
ir[N] = 1  # Right of rightmost is leftmost in a ring
il[1] = N  # Left of leftmost is rightmost

"""
Simulate the species over the given time interval.
Return the age of each species at each time step and 
the minimum fitness factor at each time step.
"""
function simulate!(species, t_space)
    age = zeros(N, t_end)
    mins = zeros(t_end)
    
    for t in t_space[2:end]
        age[:, t] = age[:, t-1] .+ 1
        # We want to find the weakest species.
        # `findmin` returns the minimum value and the index of that value.
        mins[t], index_min = findmin(species)
        species[index_min] = rand()
        species[ir[index_min]] = rand()  # Right
        species[il[index_min]] = rand()  # Left
        
        age[index_min, t] = 0
        age[ir[index_min], t] = 0  # Right
        age[il[index_min], t] = 0  # Left
    end
    
    return age, mins
end

age, mins = simulate!(species, t_space);
# We add a semicolon so that Jupter does not print the values of age and mins.

using PyPlot
# There is a lot happening in this first line. Let's break it down.
# The `[:, 1:100:end]` resembles Python's slicing syntax, but there are slight differences.
# The first colon simply means every row, just like Python.
# `1:100:end` follow the syntax `START:STEP:END`, so it means
# 'from first column to last column, take every 100th column'.
# `transpose` obviously transposes the matrix, this is done only
# to have the the axis where we think they make the most sense,
# ie. time on second axis.
# `vmax` is the same as in matplotlib. It is included to
# give better contrast in the figure, by setting all ages over
# `vmax` to the same color, bringing out the details in
# rapidly changing areas.
plt.imshow(transpose(age[:, 1:100:end]), vmax=30000)
plt.xlabel("Species")
plt.ylabel("Age [# iterations]")
plt.title("The age of species (x-axis) through the simulation.")
plt.show()

basket_size = 500
plt.plot([maximum(a) for a in Iterators.partition(mins, basket_size)])
plt.ylabel("Minimum fitness factor")
plt.xlabel("Time")
# Since we collected our results in baskets, the numbers on the time-axis are wrong.
# We either have to rescale them or remove them.
# Since we are mainly interested in the qualitative behavior, we remove the numbers.
plt.xticks([]);