Scientific computing

(for the rest of us)

Advanced iteration

In this module, we will see how we actually iterate over objects in Julia. Although the content of the previous module is very important, as it forms the basis of all ways to iterate, there are a number of functions that greatly facilitate our task. We finish this module by simulating a simple host-parasitoid model.

What are the numbers between 12 and 17? We can represent this as a UnitRange (a memory efficient way to store sequences):

our_seq = 12:17
12:17

After reading the module on indexing, we could get the second entry in this sequence with

our_seq[2]
13

But of course, performing any operation on all numbers in the sequence would be a little time-consuming. Therefore, we will iterate.

It would be very tempting to iterate from the first index (1) to the last index of the sequence (its length). Julia discourages this, because length is not really meant to help with iteration. And luckily, there is a much, much better

for i in eachindex(our_seq)
    @info our_seq[i]
end
[ Info: 12
[ Info: 13
[ Info: 14
[ Info: 15
[ Info: 16
[ Info: 17

What is this mysterious eachindex?

eachindex(our_seq)
Base.OneTo(6)

In short, it is an object that Julia prepares for us, that makes iteration safe. But there is an even better way to iterate. Let’s assume that we have a sequence of evenly spaced numbers:

our_other_seq = LinRange(0.0, 1.0, 6)
6-element LinRange{Float64, Int64}:
 0.0,0.2,0.4,0.6,0.8,1.0

We can iterate on these values to take, for example, their square, this way:

for i in eachindex(our_other_seq)
    @info "i = $(i)   (xᵢ)² = $(our_other_seq[i]^2.0)"
end
[ Info: i = 1   (xᵢ)² = 0.0
[ Info: i = 2   (xᵢ)² = 0.04000000000000001
[ Info: i = 3   (xᵢ)² = 0.16000000000000003
[ Info: i = 4   (xᵢ)² = 0.36
[ Info: i = 5   (xᵢ)² = 0.6400000000000001
[ Info: i = 6   (xᵢ)² = 1.0
Ah, yes. About this. 0.2^2.0 is not 0.04. There is a reason for this, and it is: computers are not very good with numbers. It’s OK, neither are we; hopefully it’ll sort itself out (it won’t).

But there is a more efficient way to iterate:

for (i, x) in enumerate(our_other_seq)
    @info "i = $(i)   (xᵢ)² = $(x^2.0)"
end
[ Info: i = 1   (xᵢ)² = 0.0
[ Info: i = 2   (xᵢ)² = 0.04000000000000001
[ Info: i = 3   (xᵢ)² = 0.16000000000000003
[ Info: i = 4   (xᵢ)² = 0.36
[ Info: i = 5   (xᵢ)² = 0.6400000000000001
[ Info: i = 6   (xᵢ)² = 1.0

The enumerate function is making things a little more complex because our mental model of for, variable, values is a little bit invalidated now. This is because it returns not one value but two: the position in the object we are iterating over, and the value at this position. This is a little confusing, so let’s open-up the enumerate function:

collect(enumerate(our_other_seq))
6-element Vector{Tuple{Int64, Float64}}:
 (1, 0.0)
 (2, 0.2)
 (3, 0.4)
 (4, 0.6)
 (5, 0.8)
 (6, 1.0)

This is something we know! It’s a tuple! It’s tuples in a vector! And we know from the module on tuples that they can be used to store values until we are ready to use them, and so this is what enumerate does: it stores together the position and the value.

But what about arrays with higher dimensions? The same logic applies. Let’s create a little matrix, and see how we can iterate over it:

A = reshape(Array(7:12), 2, 3)
2×3 Matrix{Int64}:
 7   9  11
 8  10  12

Let’s start to get a sense of the output of eachindex:

collect(enumerate(A))
2×3 Matrix{Tuple{Int64, Int64}}:
 (1, 7)  (3, 9)   (5, 11)
 (2, 8)  (4, 10)  (6, 12)

This is very similar to the output we got for a vector, with the exception that the shape of the enumerated elements matches the shape of the matrix. Will it be an issue? Is there something we need to do? No.

Recall from the module on indexing that we can index a matrix linearly, so we don’t need to change the way we work:

for (pos, val) in enumerate(A)
    @info "A[$(pos)] = $(val)"
end
[ Info: A[1] = 7
[ Info: A[2] = 8
[ Info: A[3] = 9
[ Info: A[4] = 10
[ Info: A[5] = 11
[ Info: A[6] = 12

But what if we wanted to use the fact that matrices have rows and columns? In this case, we can use the axes function:

axes(A)
(Base.OneTo(2), Base.OneTo(3))

When called on an array, it will return a tuple of iterators (OneTo is a weird object, but essentially, OneTo(3) will return the numbers from 1 to 3), one for each dimension. The axes function has additional methods where we specify the arguments, so we can write, for example:

for row in axes(A, 1)
    for col in axes(A, 2)
        @info "A[$(row),$(col)] = $(A[row,col])"
    end
end
[ Info: A[1,1] = 7
[ Info: A[1,2] = 9
[ Info: A[1,3] = 11
[ Info: A[2,1] = 8
[ Info: A[2,2] = 10
[ Info: A[2,3] = 12

But wait! This is two loops, one nested in the other. There has got to be an easier way to write this. When we are are dealing with nested loops, we can declare all of them on the same line:

for row in axes(A, 1), col in axes(A, 2)
    @info row, col
end
[ Info: (1, 1)
[ Info: (1, 2)
[ Info: (1, 3)
[ Info: (2, 1)
[ Info: (2, 2)
[ Info: (2, 3)