Overflow with Iterators.product

I can’t reproduce your error. Your posted code works fine for me. I also tried putting your code into a function to make it a bit faster, and accumulated a sum of the data in the loop to make sure the compiler doesn’t optimize it away:

function f(data, L, d)
    s = zero(eltype(data))
    for n in Iterators.product(repeat([1:L],d)...)
        s += data[n...]
    end
    return s
end

and this also runs without error for me on your example inputs L = 300; d = 3; data = randn((L,L,L));.

Are you sure you didn’t accidentally change L or d to be inconsistent with data?

Much, much better ways. First, for efficiency, you really want the number of loops (d) to be a compile-time constant here, and to use tuples rather than arrays so that the compiler knows the length. A simple modification of your code is to use:

Iterators.product(axes(data)...)

instead. However, it’s even cleaner to use

CartesianIndices(data)

instead — that’s what CartesianIndices are for. Note that iterating over CartesianIndices returns a CartesianIndex, not a tuple, but you can efficiently convert it to a tuple via Tuple(index). In particular, you can do:

for n in CartesianIndices(data)
    s += data[n]   # no splatting required
end

or:


for n in CartesianIndices(data)
    t = Tuple(n)
    s += data[t...]   # splatting a tuple of indices
end

Quantitatively, your original repeat method gives (via BenchmarkTools.jl):

julia> @btime f($data, $L, $d);
  6.186 s (162000023 allocations: 5.63 GiB)

due to the type instability of your runtime-dependent number of loops (d). Using Iterators.product(axes(data)...) gives me a 200× speedup:

julia> @btime f($data);
  31.091 ms (0 allocations: 0 bytes)

and using CartesianIndices(data) gives about the same performance, but more elegant code:

julia> @btime f($data);
  31.013 ms (0 allocations: 0 bytes)

PS. Please quote your code in the future to make it easier to read and to copy-paste. (I edited your post to fix this.) See PSA: how to quote code with backticks

2 Likes