Surprising dot broadcast result when indices are repeated?

General Usage
, , ,
1

I am confused about the following code. I have an indexing array pixmap to select which elements to add to:

julia> n = 10
julia> pixmap = fill(1, n)
julia> v = fill(0, n)
julia> @. v[pixmap] += 1
julia> v
10-element Vector{Int64}:
 1
 0
 0
 0
 0
 0
 0
 0
 0
 0

I find this outcome surprising. I would have expected v[1]=10. Can @__dot__ not handle indexing to the same element multiple times?

I am on julia-1.6.1 on WSL2/Ubuntu.

1 Like
2

Since you directly indexing with array you get a temporary Array for v[pixmap] in rhs so it is adding 1 to 1 10 times and stores at location 1. You need to use views to affect both sides of the assignment.

3

The reason is that:

v[pixmap] += 1

is in this case the same as:

copyto!(view(v, pixmap), v[pixmap] .+ 1)

so you have 1 set to index 1 of v 10 times.

For the same reason you have:

julia> v = fill(0, n)
10-element Vector{Int64}:
 0
 0
 0
 0
 0
 0
 0
 0
 0
 0

julia> @. v[pixmap] += 1:10
10-element view(::Vector{Int64}, [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) with eltype Int64:
 10
 10
 10
 10
 10
 10
 10
 10
 10
 10

julia> v
10-element Vector{Int64}:
 10
  0
  0
  0
  0
  0
  0
  0
  0
  0

(i.e. the last computed value is eventually stored in v)

1 Like
4

I do not think using views would resolve the issue as broadcasting will do unaliasing and the result would be the same. Could you post what you exactly meant? Thank you!

5

When I use @views @. v[pixmap] += 1. I get 10 as expected for this case but you are right there are cases where there is alias check but not for this one though.

6

Ah, you are right:

julia> n = 10;

julia> pixmap = fill(1, n);

julia> v = fill(0, n);

julia> v1 = @view v[pixmap];

julia> v2 = @view v[pixmap];

julia> @. v1 .= v2 + 1;

julia> v
10-element Vector{Int64}:
 10
  0
  0
  0
  0
  0
  0
  0
  0
  0

julia> v = fill(0, n);

julia> v1 = @view v[pixmap];

julia> v2 = @view v[copy(pixmap)];

julia> @. v1 .= v2 + 1;

julia> v
10-element Vector{Int64}:
 1
 0
 0
 0
 0
 0
 0
 0
 0
 0

and this difference is caused by:

broadcast_unalias(dest, src) = dest === src ? src : unalias(dest, src)

definition (the alias check is triggered, but in this case Julia Base decides not to perform unaliasing).

@mbauman - this is a corner case, but the question is if we want to accept it as intended, or something should be fixed here?

7

I believe this might be fixed by RFC: Fix #31392, unaliasing of broadcast arguments against destinations with repeated indices by mbauman · Pull Request #31407 · JuliaLang/julia · GitHub

8

Yes - this should fix this.

@hsgg - after this fix what @tomaklutfu proposes will not work as it does now and you will get 1 always (if I understand things correctly).

9

Thanks @bkamins @mbauman @tomaklutfu ! Writing it as (#31392)

v[[1,1]] .= [1,2]

makes it obvious that the correct answer is not obvious.

I think this is an acceptable fix, and people like me will have to write the loop explicitly.

EDIT: Made the example even simpler.

10

Thinking about this some more, the semantics I am really after is something between map() and mapreduce(), like “mapsemireduce()”. That is, consider an array of size n. Then, map() maps this array onto an array of the same size, and mapreduce() maps it onto an array of size 1, and mapsemireduce() would map it to an array of size m<n. Really, it is mapreduce() on subarrays. I don’t suppose such a programming paradigm already exists?