# Restructuring loop logic and its effect on allocations

I came across in a code base # TODO this seems very clumsy to duplicate the loop code like this? so I figured I'd give it a shot as a learning exercise. I'd like some feedback on what I came up with, particularly:

1. Is this proper idiomatic Julia, or close? Or am I just writing Python in Julia?
2. What explains the difference in allocations between the two approaches?
(Arguments used: [1,1], octaves=10)
• Old: 0.000013 seconds (63 allocations: 1.062 KiB)
• New: 0.000008 seconds (41 allocations: 1.500 KiB)

Thanks!

Original function:

# coords is [x] or [x, y] or [x, y, z] or [x, y, z, w]
function _octaves(coords::Array{T, 1} where T <: Real;
octaves::Int = 1,
persistence=1.0)
total     = 0.0
frequency = 1.0
amplitude = 1.0
maxval    = 0.0
# TODO this seems very clumsy to duplicate the loop code like this?
l = length(coords)
if l == 1
for i in 1:octaves
total += simplexnoise(coords[1] * frequency) * amplitude
maxval += amplitude
amplitude *= persistence
frequency *= 2
end
elseif l == 2
for i in 1:octaves
total += simplexnoise(coords[1] * frequency, coords[2] * frequency) * amplitude
maxval += amplitude
amplitude *= persistence
frequency *= 2
end
elseif l == 3
for i in 1:octaves
total += simplexnoise(coords[1] * frequency, coords[2] * frequency, coords[3] * frequency) * amplitude
maxval += amplitude
amplitude *= persistence
frequency *= 2
end
elseif l == 4
for i in 1:octaves
total += simplexnoise(coords[1] * frequency, coords[2] * frequency, coords[3] * frequency, coords[4] * frequency) * amplitude
maxval += amplitude
amplitude *= persistence
frequency *= 2
end
end
end


And my refactor:

function _octaves(coords::Array{T, 1} where T <: Real;
octaves::Int = 1,
persistence=1.0)
total     = 0.0
frequency = 1.0
amplitude = 1.0
maxval    = 0.0

for _ in 1:octaves
inputs = coords .* frequency
total += simplexnoise(inputs...) * amplitude

maxval += amplitude
amplitude *= persistence
frequency *= 2
end

end


This is very idomatic, in principle. The loop refactoring is exactly what I would have done.

Here's some other suggestions:

const Coords{T} = NTuple{N, T} where N

function _octaves(coords::Coords{T}; octaves::Int=1, persistence::T=one(T)) where {T<:Real}
if octaves < 1
# prevent zero division
throw(DomainError(octaves, "octaves must be at least 1!"))
end

total = zero(T)
maxval = zero(T)

for octave in 1:octaves
p = octave - 1
frequency = 2 ^ p
amplitude = persistence ^ p

maxval += amplitude

scaled_coords = coords .* frequency
total += simplexnoise(scaled_coords...) * amplitude
end

end

1. Using one and zero to ensure type stability (a very important concept, which you should read about if you don't know it!)
2. Replacing some intermediate updates by explicit powers. The calculation of frequency is trivial to the compiler, since it is an integer power of 2. amplitude might do some unnecessary work; I think it is clearer that way, and not worth the microoptimisation, but you'd have to test that for your use case.

The choice of type for coords is really determined by the rest of the code. If this is only ever going to be used for coordinates of bounded dimension, I'd very much prefer a tuple as I wrote above. If you use vectors all the way down (this looks a bit like signal processing?), then staying with Coords{T} = Vector{T} is easier.

While not part of the question, it might be nicer to have simplenoise work on an iterable, not a varargs argument. Possibly even as simplenoise(frequency, coords), then the allocation of scaled_coords is removed.

If there's any chance that the summation will overflow, take a look at Base.widen. Also, you might want to ensure that octaves < 8 * sizeof(Int) or something like that, and add this to the check at the beginning.

If this is really the hot part of some code to be optimized, with octaves being constant over many calls, there's the possibility to write a @generated function to unroll the loop completely, but unless it is critical, it's not worth it.

As for the difference in benchmarking: first, you should do that with BenchmarkTools.jl, not @time. Then, I haven't measured myself, but it's quite surely due to the intermediate inputs array you construct.

• This is great, thank you! Looks like I have some reading lined up for tonight, I really appreciate you taking the time. – Oso May 27 '20 at 18:41