I have some code that accumulates a list of indices in a vector that need to be incremented, and then does the incrementing after the fact with vector[indices] += increment (I do it this way instead of incrementing as I accumulate because the increment depends on the number of indices being incremented). The problem arises when there are duplicate indices. I expected this code:
x = zeros(3)
indices = [2, 2]
x[indices] .+= 1
to be equivalent to this:
x = zeros(3)
indices = [2, 2]
for index in indices
x[index] += 1
end
However, the first produces [0.0, 1.0, 0.0] while the second produces [0.0, 2.0, 0.0]. I understand broadcast is supposed to be equivalent to ., so I also tried
broadcast!(v -> v + 1, x[indices], x[indices])
which does not modify x at all, though this is not so surprising since x[indices] is creating a copy. Using a view reproduces the behavior of the . operator:
broadcast!(v -> v + 1, @view(x[indices]), x[indices])
It looks like this behavior is consistent all the way back to Julia 1.0. Is this the expected behavior?
If you’re interested in the actual code that led me to discover this, see here
There are two things going on here. One is that lhs .+= 1 is the very simple transform that expands to lhs .= lhs .+ 1. And what you have on the left hand side is x[indices] — which, when moved onto the right-hand side, becomes the standard copying indexing behavior. So you end up with:
x[[2,2]] .= x[[2,2]] .+ 1
Now if you’re enterprising, you can remember that @views will avoid making that copy! So you can do @views x[[indices]] .+= 1 — and indeed that returns [0.0, 2.0, 0.0]. But now you’re at the edge of hitting the second problem, and that’s that Julia’s semantics here simply aren’t very well defined. Broadcast makes a “best-effort” attempt at avoiding strange aliasing behaviors, but it’s not 100% clear what the behaviors should be. This is How should in-place broadcasting with repeated indices behave? · Issue #31392 · JuliaLang/julia · GitHub.
1 Like
Makes sense, thanks! FWIW, numpy has the same behavior.
The part that makes it particularly quirky (and likely to be troublesome) is that the unaliasing behaviors depend upon the exact identity of your indices:
julia> indices = [1,1];
julia> x = [3]; @views x[indices] .= x[[1,1]] .+ x[[1,1]]; x
1-element Vector{Int64}:
6
julia> x = [3]; @views x[indices] .= x[indices] .+ x[[1,1]]; x
1-element Vector{Int64}:
9
julia> x = [3]; @views x[indices] .= x[indices] .+ x[indices]; x
1-element Vector{Int64}:
12
julia> x = [3]; @views x[[1,1]] .= x[indices] .+ x[indices]; x
1-element Vector{Int64}:
6
I wish we would’ve identified this before 1.0; then I would’ve pushed harder to fix it. But now fixes are much tricker because some folks expect (as you did) that it’s an iterative method while others think of it as a whole-array-at-once method.
1 Like
I wonder if a useful “fix” would be to just issue a warning if there are duplicate indices (maybe that’s difficult/computationally expensive to figure out though). Kinda like the pandas SettingWithCopyWarning, then at least I would have checked this and figured it out earlier.