Should a = [0]; a[[1,1]]+=[1,2] cause a==[3]?

Internals & Design
,
1

To repeat the title, should a = [0]; a[[1,1]]+=[1,2] update the array a such that a==[3]?

This is not really a pressing issue, but changing to this kind of interpretation would be helpful in some cases when dealing with certain computations on multidimensional arrays, and the current result a==[2] is somewhat biased.

I am mentioning this since manually dealing with repeating values as above in a multidimensional case can be quite bothersome, however, this might be easy to change at the fundamental level of Julia. Unfortunately, I am not able to judge about the latter part, so maybe someone can help out with that.

2

I think this is just because a[[1,1]] creates a copy.

This works:

julia> a = [0]
1-element Vector{Int64}:
 0

julia> @views a[[1,1]] .+= [1,2]
2-element view(::Vector{Int64}, [1, 1]) with eltype Int64:
 3
 3

julia> a
1-element Vector{Int64}:
 3
3 Likes
3

This is nice to know. I was not aware of the macro, thank you.
Though should the command not always give a view when used as left value?

I find some of the following behavior still slightly unexpected. Maybe you know the underlying reasoning of the following commands.

Technically, you can probably explain why the current Julia (or rather 1.7.2) behaves like this, but should it?

Macro without pointwise operation:

julia> a = [0]; @views a[[1,1]] += [1,2]
1-element Vector{Int64}:
 2

Macro with pointwise operation:

julia> a = [0]; @views a[[1,1]] .+= [1,2]
1-element Vector{Int64}:
 3

Direct view (with or without pointwise does not matter):

a = [0]; b = view(a,[1,1]); b[:] += [1,2]; a
1-element Vector{Int64}:
 2

No explicit view (with or without pointwise does not matter):

a = [0]; a[[1,1]] = a[[1,1]] .+ [1,2]
1-element Vector{Int64}:
 2

I mean, a[[1,2]] += [1,2] also changes a = [0,0] to a = [1,2].

4

This is a bit tricky on multiple levels. To tease it apart, take it one step at a time:

  • a[[1,1]] = [1,2] does an iterated indexed assignment into a. Thus it will assign 1 into a[1] and then 2 into a[1].
  • a[[1,1]] += [1,2] lowers to a[[1,1]] = a[[1,1]] + [1,2] — that is, exactly the expression on the left hand side gets placed onto the right hand side. And non-scalar indexing makes a copy by default.
  • @views a[[1,1]] += [1,2]: simply adding @views still results in a == [2] because the + operation makes a copy before doing the assignment. In other words, it doesn’t really matter whether the a[[1,1]] on the right hand side is a view or not because adding a[[1,1]] + [1,2] will happen before doing a single assignment back into a.
  • a[[1,1]] .= a[[1,1]] .+ [1,2] now we’re using .-broadcasting to “fuse” the assignment and the + operation together, but we’re still “greedily” evaluating a[[1,1]] to a new copy before doing any of it.
  • a[[1,1]] .= @views(a[[1,1,]]) .+ [1,2] will finally do the trick — and that’s functionally the same as the more concise @views a[[1,1]] .+= [1,2]. This works because now we’re finally not making any copies anywhere along the entire chain.. See below

One way to help keep this all straight is to remember that Julia will always evaluate all the arguments to a function before calling the function. There’s lots of fancy “sugar” going on here between indexed assignment (which becomes setindex! function calls), the += which just literally throws the LHS onto the RHS, and fused broadcast-., which is a beast of its own. You can also ask Julia for its intermediate work, but it takes some understanding:

julia> Meta.@lower a[[1,1]] += [1,2]
:($(Expr(:thunk, CodeInfo(
    @ none within `top-level scope`
1 ─ %1 = Base.vect(1, 1).          # construct the literal [1,1]
│   %2 = Base.getindex(a, %1).     # do a[[1,1]] — making a copy
│   %3 = Base.vect(1, 2)           # construct the literal [1,2]
│   %4 = %2 + %3                   # add the two vectors, creating a copy
│        Base.setindex!(a, %4, %1) # do the indexed assignment
└──      return %4
))))
7 Likes
5

Okay, I am getting there. I am a bit used to C++ style where += may rather not create any copies.

But testing a[[1,1]] .= @views(a[[1,1,]] .+ [1,2]), it does not work for me. It yields [2] as does @views a[[1,1]] .= a[[1,1]] .+ [1,2].

Really the only thing that seems to work is @views a[[1,1]] .+= [1,2] which does indeed yield [3].

1 Like
6

Ooof, yeah, you’re right. I think we’re now running afoul of Broadcast’s aliasing detection which tries to ensure we don’t end up with odd situations like this. This is indeed an odd case — and worth an issue.

This goes well into the weeds now, but the issue is that [1,1] !== [1,1] — each new invocation of [1,1] is a newly created vector… and due to some buggy heuristics this matters for the aliasing detection.

julia> a = [0]; @views a[[1,1]] .= a[[1,1]] .+ [1,2]; a
1-element Vector{Int64}:
 2

julia> a = [0]; idx = [1,1]; @views a[idx] .= a[idx] .+ [1,2]; a
1-element Vector{Int64}:
 3
3 Likes
7

I think the following captures it quite well:

Okay, firstly, the following works

a = [0]
b = view(a,[1,1]) # or Base.dotview(b,[1,1])
b .= b .+ [1,2]

in the sense that we now have a == [3]. Naturally, just b = ... would just assign to b whatever is right of the equal sign. That makes sense.

However, this is where I get seriously confused and I think this is what you pointed out. Looking at the underlying procedure, we have

julia> Meta.@lower b .= b .+ [1,2]
:($(Expr(:thunk, CodeInfo(
    @ none within `top-level scope`
1 ─ %1 = Base.vect(1, 2)
│   %2 = Base.broadcasted(+, b, %1)
│   %3 = Base.materialize!(b, %2)
└──      return %3
))))

Here %3 = Base.materialize!(b, %2) does some magic broadcasting the assignment operation together with the addition. I guess because it first adds, then assigns, then adds, then assigns, instead of twice adding, then twice assigning.

Buuut… now take this:

julia> Meta.@lower a[[1,1]] .= b .+ [1,2]
:($(Expr(:thunk, CodeInfo(
    @ none within `top-level scope`
1 ─ %1 = Base.vect(1, 1)
│   %2 = Base.dotview(a, %1)
│   %3 = Base.vect(1, 2)
│   %4 = Base.broadcasted(+, b, %3)
│   %5 = Base.materialize!(%2, %4)
└──      return %5
))))

The first two steps are just what we did earlier for b. The last three steps are literally the three steps above, however, we then end up with a == [2]? This is seriously weird, if I did not just miss something.
I get that we now have two different views a[[1,1]] (as appearing as L-value) and b, so a!===b, but is this how things should work out?

What surprisingly works as well is view(a,[1,1]) .+= [1,1] (I made an error here ealier). We get a==[3].
But what does not work is view(a,[1,1]) .= view(a,[1,1]) .+ [1,1]. Here we get a==[2], because the first and the second [1,1] are not === equal. Yeah, it can not be right :sweat_smile:. As it seems, for the first call, [1,1] is actually evaluated first, and then += is interpreted and expanded :face_with_raised_eyebrow:.

The rest of this post is drifting away from the initial discourse a bit, but maybe I can explain my whole confusion with += and .+=

While I have no idea how to deal with the above, the literal expansion of x+=y really makes my head hurt even though technically it is simple. Take this:

a = [0]
b = view(a,[1,1])
b += [1,2]

Here, a==[0], as after the expansion to b = b + 1, the left b is now just assigned as the new vector [1,2] because it is no longer acting as view. To do so, one has to use b .+= [1,2], which however at the same times changes the order of assigning and adding, and we get a = [3].

On the other hand, with

a = [0]
a[[1,1]] .+= [1,2]

we now have (really just the other way around) a view on the left, but when expanded to the right, now a[[1,1]] is no longer a view, and as carried out above, we get a==[2]. This kind of makes me want to avoid += altogether.

Is there even a particular reason why a[[1,1]] as well as a[:] is not just always a view? I mean, the adjoint a' is a view, reshapings are views, but vectorization a[:] and subarrays are not…? If so, then a[[1,1]] .+= [1,2] (modulo aboves weird behavior) would indeed yield [3] even without further change of +=.

1 Like
8

(I am assuming you meant to write [1,1] !=== [1,1]):

The operators are:

== and !=
=== and !==
9

(Ah, sorry. Did not even think about it. Of course.)

1 Like