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