`broadcast` for `AbstractArray` types

[ettersi]

Consider the following code:

immutable MyArray{T,N} <: AbstractArray{T,N}
    a::Array{T,N}
end
Base.linearindexing{T,N}(::Type{MyArray{T,N}}) = Base.LinearFast()
Base.size(m::MyArray) = size(m.a)
Base.getindex(m::MyArray, i::Int) = m.a[i]
Base.setindex!(m::MyArray, v, i::Int) = m.a[i] = v
Base.similar{T,N}(::MyArray, ::Type{T}, dims::NTuple{N,Int}) = MyArray(Array(T,dims))

Thanks to the AbstractArray magic, these 8 lines should define a type with the same interface as Array. This is however not the case:

julia> typeof(MyArray([1]).*MyArray([1]))
Array{Int64,1}

The problem originates from the following definition of broadcast():

broadcast(f, As...) = 
    broadcast!(
        f, 
        Array(promote_eltype(As...), broadcast_shape(As...)),
        As...
    )

I suggest changing this definition to:

broadcast(f, As...) = 
    broadcast!(
        f, 
        similar(first(As), promote_eltype(As...), broadcast_shape(As...)),
        As...
    )

This would fix the above problem.

[pabloferz]

#17389 should fix this.

[yuyichao]

No it doesn't?

WARNING: Error during initialization of module Random:
Base.MethodError(f=Base.#promote_array_type(), args=(Base.#.+(), Array{UInt32, 1}[0x640c92bc, 0x07683af9, 0xe5e6e6d2, 0xb346b00e], 0x00000001))
               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.5.0-dev+12987 (2016-07-21 12:37 UTC)
 _/ |\__'_|_|_|\__'_|  |  pz/pobc/70a6696 (fork: 2 commits, 0 days)
|__/                   |  x86_64-pc-linux-gnu

julia> immutable MyArray{T,N} <: AbstractArray{T,N}
           a::Array{T,N}
       end

julia> Base.linearindexing{T,N}(::Type{MyArray{T,N}}) = Base.LinearFast()

julia> Base.size(m::MyArray) = size(m.a)

julia> Base.getindex(m::MyArray, i::Int) = m.a[i]

julia> Base.setindex!(m::MyArray, v, i::Int) = m.a[i] = v

julia> Base.similar{T,N}(::MyArray, ::Type{T}, dims::NTuple{N,Int}) = MyArray(Array(T,dims))

julia> typeof(MyArray([1]).*MyArray([1]))
Array{Int64,1}

[pabloferz]

@yuyichao You're right. I misred the OP. And I agree that this should use similar instead of Array.

[stevengj]

similar is guaranteed to be mutable, and seems like the right choice here.

[tkelman]

Don't we currently have map always returning an Array, and map being inconsistent with broadcast might be strange, wouldn't it?

[stevengj]

Good point, @tkelman. Note that if you want to control the output type, you can always pre-allocate the output and call broadcast! (or use .= after #17510).

[timholy]

I think it's inconsistent to make the output type depend on only the first of the arrays. Might need a typejoin---but that could introduce non-leaftypes.

[ettersi]

I agree with what timholy says, but in fact it's what Julia does for binary operations. Try

julia> typeof(MyArray([1])+[1])
MyArray{Int64,1}

julia> typeof([1]+MyArray([1]))
Array{Int64,1}

To overcome this, you would probably have to introduce a new "promotion" system for AbstractArrays. Not sure whether that's worth the effort.

[ssfrr]

@tkelman map doesn't currently return an Array, I think it uses similar:

julia> a = MyArray(rand(4, 2))
4×2 MyArray{Float64,2}:
 0.981792  0.116486
 0.232802  0.519291
 0.736532  0.234633
 0.546964  0.80054

julia> map(x->x*2, a)
4×2 MyArray{Float64,2}:
 1.96358   0.232972
 0.465605  1.03858
 1.47306   0.469265
 1.09393   1.60108

I agree that basing the output type on the first arg seems weird. I've been thinking about a promote_similar method, but I'm not sure how you'd define it to dispatch correctly on a Tuple that contains at least one MyArray...

Then we could do something like:

broadcast(f, As...) =
    broadcast!(
        f, 
        promote_similar(As, promote_eltype(As...), broadcast_shape(As...)),
        As...
    )

Possibly promote_similar could take the type as the first arg, like:

promote_similar(::Type{MyArray}, As, T, dims)

and then broadcast would be:

function broadcast(f, As...)
    T = promote_type(map(typeof, As)...)
    broadcast!(
        f, 
        promote_similar(T, As, promote_eltype(As...), broadcast_shape(As...)),
        As...
    )
end

Feels pretty fiddly. :/

I'm definitely interested in finding a solution though. I hit this trying to make sure I can do arithmetic on my SampleBuf types in SampledSignals.jl.

[pabloferz]

With the current state of affairs, we have now

(MyArray([1]) + [1])::Vector{1}
([1] + MyArray([1]))::Vector{1}

that is, you now get an Array. If you want a MyArray you have to write Base.Broadcast.promote_containertype methods for your own type. You can find an example of this in the tests and a probable simpler future version on #20102.

I believe that should solve this.

[stevengj]

If you want a specific container type for the output, you can also just use .= (i.e. broadcast!)

[tkelman]

Given the number of operations that go through broadcast now, I don't think requiring in place operations to get these to return a specific array type is a good enough answer for extensibility.

[martinholters]

Also, maybe one wants to use the broadcast machinery for determining the element type, and hence cannot pre-allocate the destination container.

[mbauman]

I've been thinking about this issue with respect to the broadcast revamp. Things are much improved since its inception given that we now have a documented interface for extending broadcast. It's still somewhat unsatisfactory, though, in that we don't yet have map(f, A) == f.(A), and customizing the result type there is relatively complex.

One possible solution would be for the fallback DefaultArrayStyle implementation to use similar(A, ...) in cases where there is only one array argument. Of course, the very next thing you'll want is A .* 2 or perhaps even A .+ A. I neither see a clear line for delineating such a special case nor a general rule that gracefully collapses to the one-argument behavior we want.