`map!` does not check for aliasing

[jakobnissen]

julia> x = [1,2]
2-element Vector{Int64}:
 1
 2

julia> map!(i -> i + 1, x, view(x, 2:-1:1))
2-element Vector{Int64}:
 3
 4

This is silently incorrect. The result should be [3, 2].
Throwing an error would probably be preferable.

[roflmaostc]

This seems to work however:

julia> x = [1,2]
2-element Vector{Int64}:
 1
 2

julia> y = view(x, 2:-1:1)
2-element view(::Vector{Int64}, 2:-1:1) with eltype Int64:
 2
 1

julia> map!(i  -> i+1, y, y)
2-element view(::Vector{Int64}, 2:-1:1) with eltype Int64:
 3
 2

julia> y
2-element view(::Vector{Int64}, 2:-1:1) with eltype Int64:
 3
 2

julia> x
2-element Vector{Int64}:
 2
 3

[MasonProtter]

Broadcast has the same problem:
Broadcast does not have this problem:

julia> x = [1,2]
2-element Vector{Int64}:
 1
 2

julia> x .= (i -> i + 1).(view(x, 2:-1:1))
2-element Vector{Int64}:
 3
 2

julia> x
2-element Vector{Int64}:
 3
 2

[roflmaostc]

Broadcast has the same problem:

julia> x = [1,2]
2-element Vector{Int64}:
 1
 2

julia> x .= (i -> i + 1).(view(x, 2:-1:1))
2-element Vector{Int64}:
 3
 2

julia> x
2-element Vector{Int64}:
 3
 2

Isn't that example ok?
I mean, we want to increment the reverse array of x and assign it to x (which is a vector)?
We can also see that there is an allocation probably to avoid the aliasing?

julia> @time x = collect(1:1000);

julia> f = (i -> i + 1);

julia> @time x .= f.(view(x, 1000:-1:1));
  0.054761 seconds (123.19 k allocations: 6.564 MiB, 15.81% gc time, 99.45% compilation time)

julia> @time x .= f.(view(x, 1000:-1:1));
  0.000027 seconds (4 allocations: 8.125 KiB)

[MasonProtter]

Isn't that example ok?

🤦 yeah, I meant to write the opposite. Broadcasting does not have this problem. I've edited my comment to fix that, thanks

[CameronBieganek]

mul! has the following warning in its docstring: "Note that Y must not be aliased with either A or B." Maybe the map! docstring needs a similar warning.

Throwing an error would probably be preferable.

Off the top of my head, I'm not even sure if it's possible to detect arbitrary kinds of aliasing, other than "egal": y === x. Actually, a view probably does't count as an alias. An alias is just another name for the same thing. So, I'm coming around to the opinion that the behavior demonstrated by OP is correct, even though it is counter-intuitive.

Note that map! works fine for actual aliasing:

julia> y = x = [4.0, 9.0];

julia> map!(sqrt, y, x)
2-element Vector{Float64}:
 2.0
 3.0

Perhaps the docstring for map! can have some clarification along these lines:

The transformation of collection and the storage into destination happens sequentially. If collection and destination are dependent in some way, then the results might be counter-intuitive.

[roflmaostc]

Why is it counter intuitive?

Look at the difference closely, it's 2 and not 1. So not only the order might be wrong but also the resulting numbers are wrong.

[jakobnissen]

@CameronBieganek "Aliasing" here refers to memory aliasing, i.e. having two variables which point to the same memory.

The behaviour in OP is absolutely not correct. The number 4 is in the output array, and 4 is neither 1+1 nor 2+1, so the behaviour as described in the documentation does not hold. The underlying implementation, namely a single loop over the elements which both reads and writes, leaks through.

I'm not sure clarifying that the transformation happens sequentially is really a help here - surely, that is just an implementation detail. For example, if the loop SIMDs or is parallelized, the sequential behaviour may be violated, and we would not want to disallow that.

The best solution would be to somehow statically detect aliasing and throw an error. I'm not a compiler guy - I have no idea if that is feasible. I suspect it's not, without some Rust-like borrowchecker. If it is not possible to detect aliasing in a way which is both reliable an fast, then having a warning in the docs is the next-best idea.

Julia does provide Base.mightalias, but it just returns a pointer to the start of the array (or parent array, in case of views), so it's quite conservative. Throwing an error if you try to map! two aliasing arrays would be breaking.

[CameronBieganek]

Look at the difference closely, it's 2 and not 1.

A little respect, please. I'm not blind. If the intended semantics of map! are that f is applied sequentially, then the output is correct.

so the behaviour as described in the documentation does not hold.

This is the problem. The docstring for map! is so vague that it's impossible to know which of the two possible behaviors is intended. Basically we need an authority figure to chime in and tell us which of the following two semantic definitions for map! is the intended one:

1. First transform collection by applying f to each element, and then store the result in destination.
2. Sequentially apply f to collection, with results sequentially stored in destination.

Given the possibility of memory aliasing and given that map! is supposed to operate in-place, I think definition 2) is the only reasonable one.

Note that the docstring for map can get away with being a little vague, since there's no mutation involved, but for map! we need to be more precise.

[jakobnissen]

In my view it's clear that the loop order and load/store order is an implementation detail. For example, if, in the OP example, two values were loaded at a time, then both were stored it would not error. This kind of loop unrolling is a completely standard transformation.

It would be highly strange if load/store order was part of the documented API. That seem like such a low-level detail and seems odd to expose at the API level. Further, if it is exposed, then the current implementation is still buggy, since we must guarantee that it never SIMDs.

[CameronBieganek]

It would be highly strange if load/store order was part of the documented API. That seem like such a low-level detail and seems odd to expose at the API level.

There's nothing odd about exposing sequential behavior at the API level. Look at foldl, foldr, and just about anything related to iteration. I think "sequentially" and "iteratively" are close to synonymous, and there's plenty of iteration in the language. How about I change the wording of definition 2) to the following?

2. Iteratively apply f to collection. At each iteration, store the return value of f at the corresponding location in destination.

As I mentioned already, due to the mutation involved, the semantics of map! needs to be made more precise. And due to the in-place nature of map!, the semantics of definition 1) are not really tenable.

Further, if it is exposed, then the current implementation is still buggy, since we must guarantee that it never SIMDs.

As far as I'm aware, there is no interface for declaring that SIMD is not allowed, so we have to rely on unit tests. Here's the code for the primary map! method:

function map!(f::F, dest::AbstractArray, A::AbstractArray) where F
    for (i,j) in zip(eachindex(dest),eachindex(A))
        val = f(@inbounds A[j])
        @inbounds dest[i] = val
    end
    return dest
end

I'm not a compiler guy either, but the docstring for @simd says "The object iterated over in a @simd for loop should be a one-dimensional range." Of course, the loop in map! is not annotated with @simd, but this leads me to believe that the Julia compiler would not attempt SIMD on this loop, since the object being iterated over is an Iterators.Zip object, not a range object. Even if the object being iterated over was a range, I wonder if this is one of the situations where the compiler recognizes that it is not safe to change the order of iteration, and therefore does not attempt SIMD. Anyhow, with a decent test suite for this behavior (which we don't have) and CI/CD on multiple platforms (which we do have), we can be reasonably confident that SIMD is not happening.

[LilithHafner]

I think we should continue to support egal map!(f, x, x) and other non-offset aliasing map!(f, x, view(x, 1:3)), and of course non-aliased map!(f, x, [1,2,3]). When there is offset aliasing (e.g. the OP), I think we should declare undefined behavior to support simd, parallelization, etc. How does this docstring look?

---

`map!(function, destination, collection...)`

Like map, but stores the result in destination rather than a new
collection.

destination must be at least as large as the smallest collection.

If destination shares memory with but is not egal to (===) any of the
collections, behavior may be undefined. See extended help for more info.

See also: map, foreach, zip, copyto!.

Examples

julia> a = zeros(3);

julia> map!(x -> x * 2, a, [1, 2, 3]);

julia> a
3-element Vector{Float64}:
 2.0
 4.0
 6.0

julia> map!(+, zeros(Int, 5), 100:999, 1:3)
5-element Vector{$(Int)}:
 101
 103
 105
   0
   0

Extended Help

The load-store order used by map! is an implementation detail. For example, it may
sequentially load from the collections and store to destination in one element at a time
or it may load multiple elements from the collections, perform computation in parallel and
then store multiple elements to destination.

Because of this, if the contents of the ith index of destination impacts the value
of the jth index of a collection where i ≠ j and j ∈ eachindex(destination), the
behavior of map! is undefined.

Examples

julia> a = [1:4;];

julia> map!(x -> x+1, a, a) # no offset, so okay
4-element Vector{Int64}:
 2
 3
 4
 5

julia> map!(x -> x+10, a, view(a, 3:4)) # no overlap, so okay
4-element Vector{Int64}:
14
15
 4
 5

julia> map!(x -> x+100, a, view(a, 1:3)) # no offset, so okay
3-element Vector{Float64}:
114
115
104
  5

julia> map!(x -> x+100, a, view(a, 2:4)) # undefined
[...]

julia> map!(x -> x+100, view(a, 2:4), view(a, 1:3)) # undefined
[...]

julia> map!(x -> x+100, a, view(a, 4:-1:1)) # undefined
[...]