Restructuring loop logic and its effect on allocations

Question

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
    return total / maxval
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

    return total / maxval
end

Answer 1

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

    return total / maxval
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.