fix #4883, result type of broadcast for arbitrary functions

[JeffBezanson]

Silly demo:

julia> broadcast(tuple, [1 2 3], ["a", "b", "c"])
3×3 Array{Tuple{Int64,String},2}:
 (1,"a")  (2,"a")  (3,"a")
 (1,"b")  (2,"b")  (3,"b")
 (1,"c")  (2,"c")  (3,"c")

julia> broadcast((x,y)->(rand(Bool)?1.0:x,y), [1 2 3], ["a", "b", "c"])
3×3 Array{Tuple{Real,String},2}:
 (1.0,"a")  (2,"a")    (1.0,"a")
 (1,"b")    (1.0,"b")  (1.0,"b")
 (1.0,"c")  (1.0,"c")  (1.0,"c")

This still uses promote_eltype_op, but falls back to the "map" algorithm when that returns Any.

Reducing the size of this code is a good challenge for the future. The problems are performance, and handling index ranges and iteration order of the output properly in collect. Another issue is that you can't use @simd on a while loop.

[stevengj]

Still doesn't seem fully fixed:

julia> broadcast(sin, [1,2,3])
ERROR: InexactError()
 in macro expansion at ./broadcast.jl:97 [inlined]
 in macro expansion at ./simdloop.jl:73 [inlined]
 in macro expansion at ./broadcast.jl:91 [inlined]
 in _broadcast!(::Base.#sin, ::Array{Int64,1}, ::Tuple{Tuple{Bool}}, ::Tuple{Array{Int64,1}}, ::Type{Val{1}}) at ./broadcast.jl:86
 in broadcast! at ./broadcast.jl:139 [inlined]
 in broadcast_t at ./broadcast.jl:195 [inlined]
 in broadcast(::Function, ::Array{Int64,1}) at ./broadcast.jl:197
 in eval(::Module, ::Any) at ./boot.jl:234
 in macro expansion at ./REPL.jl:92 [inlined]
 in (::Base.REPL.##1#2{Base.REPL.REPLBackend})() at ./event.jl:46

julia> map(sin, [1,2,3])
3-element Array{Float64,1}:
 0.841471
 0.909297
 0.14112 

[nalimilan]

That example works with #16995.

[stevengj]

Great! Could you add a broadcast(sin, [1,2,3]) test with a cross-reference to #4883 so that #4883 gets closed (again) when it is merged?

[stevengj]

Actually, it sounds like #16995. probably doesn't fully fix #4883, since you only fixed promote_op for numbers. For example, try length.(["foo", "bar"]).

[stevengj]

My understanding is that the fix for this shouldn't rely on promote_op ... there should be no place where map gets the correct type but broadcast does not.

[JeffBezanson]

Yes, promote_op has definitions that just guess, or that get the wrong answer, and all those should be removed. We only need a small number of promote_op definitions to keep things like + and * fast, since otherwise we lose some optimizations in that case.

[nalimilan]

@stevengj #16995 also fixes your second example. I've just added it to the tests.

[nalimilan]

@JeffBezanson Why do we lose optimizations without promote_op methods? When inference guesses a concrete type for the result, isn't that enough? Could return type declarations help?

[stevengj]

Doesn't this mean that f is called twice on the first element of the array? Can we avoid this, e.g. by passing the first element as an extra argument to _broadcast!?

Calling f twice is unexpected if you broadcast a function with side effects. Also, it's undesirable in cases where f is extremely expensive (e.g. f solves a huge PDE).

[JeffBezanson]

No; the state from next is passed to _broadcast! so it can pick up where we left off.

[JeffBezanson]

@nalimilan The reason is related to @stevengj 's question above --- we handle the first element outside the main broadcasting loop, and need a reference to the iteration state. Therefore I used a while loop, which @simd doesn't support. That's the only problem as far as I know. It's possible that (1) we could insert the simd metadata manually, (2) @simd could be easily modified to handle this, or (3) the loop is getting vectorized anyway.

[stevengj]

@JeffBezanson, with the current master broadcast calls the function twice for the first element:

julia> broadcast(println, 1:5)
1
1
2
3
4
5
5-element Array{Void,1}:
 nothing
 nothing
 nothing
 nothing
 nothing

[JeffBezanson]

Unlucky choice of arguments :)

julia> broadcast(println, 2:5)
1
2
3
4
5
4-element Array{Void,1}:
 nothing
 nothing
 nothing
 nothing

The first 1 is due to a call from promote_op.

[nalimilan]

Another consequence of that fallback: #17300 (comment)