Subtracting 0-dimensional AbstractDiskArrays fails

[sethaxen]

The problem

Here's an MWE:

julia> using DiskArrays.TestTypes

julia> v = UnchunkedDiskArray(rand(2));

julia> v-v;  # works

julia> a = UnchunkedDiskArray(fill(rand()));

julia> a - a; # fails
ERROR: MethodError: no method matching similar(::Base.Broadcast.Broadcasted{…}, ::Type{…}, ::Tuple{})
The function `similar` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  similar(::Base.ReshapedArray, ::Type, ::NTuple{N, Int64} where N)
   @ Base reshapedarray.jl:252
  similar(::Array, ::Type, ::NTuple{N, Int64}) where N
   @ Base array.jl:377
  similar(::Base.ReinterpretArray, ::Type, ::NTuple{N, Int64} where N)
   @ Base reinterpretarray.jl:223
  ...

Stacktrace:
 [1] similar(bc::Base.Broadcast.Broadcasted{DiskArrays.ChunkStyle{…}, Nothing, typeof(-), Tuple{…}}, ::Type{DiskArrays.BroadcastDiskArray{…}})
   @ Base.Broadcast ./broadcast.jl:227
 [2] broadcast_preserving_zero_d
   @ ./broadcast.jl:884 [inlined]
 [3] -(A::UnchunkedDiskArray{Float64, 0, Array{Float64, 0}}, B::UnchunkedDiskArray{Float64, 0, Array{Float64, 0}})
   @ Base ./arraymath.jl:8
 [4] top-level scope
   @ REPL[8]:1
Some type information was truncated. Use `show(err)` to see complete types.

Same error happens for ChunkedDiskArray.

Proposed solution

This happens because subtracting AbstractArrays ultimately hits this method, which for 0-dimensional arrays calls similar(bc, ::Type{T}) on the result of Base.broadcasted(-, a, a). For a normal Array, we hit these two methods, but DiskArrays lacks these analogous methods:

Base.similar(::Base.Broadcast.Broadcasted{DiskArrays.ChunkStyle{N}}, ::Type{ElType}, dims) where {N,ElType} =
    similar(Array{ElType, length(dims)}, dims)
Base.similar(::Base.Broadcast.Broadcasted{DiskArrays.ChunkStyle{N}}, ::Type{Bool}, dims) where N =
    similar(BitArray, dims)

This isn't quite enough to make the above example work, because DiskArrays also lacks an analog to this method:

@inline Base.copy(bc::Base.Broadcast.Broadcasted{DiskArrays.ChunkStyle{0}}) = bc[CartesianIndex()]

With these 3 overloads, the example works:

julia> v - v
2-element DiskArrays.BroadcastDiskArray{Float64, 1, Base.Broadcast.Broadcasted{DiskArrays.ChunkStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(-), Tuple{UnchunkedDiskArray{Float64, 1, Vector{Float64}}, UnchunkedDiskArray{Float64, 1, Vector{Float64}}}}}

Chunked: (
    [2]
)


julia> a - a
0-dimensional Array{Float64, 0}:
0.0

Does it make sense to add these?

[rafaqz]

Good catch and simple solution, seems mergeable to me.

[felixcremer]

But doesn't that change the container type of the output?
I would have expected, that a - a returns a UnchunkedDiskArray and not an Array.

[rafaqz]

a - a would normally return a BroadcastDiskArray (so the container type essentially always changes). But probably we dont need to be lazy for the zero dim case?

[felixcremer]

No we don't need to be lazy for the 0-dim case, but I am wondering about two things:

First, with this change the resulting type of operations depend on the dimensionality. I think that could be confusing.
Secondly the definition for similar on a ChunkStyle would be hit for all dimenisions and I can't oversee what this would mean for multidimensional arrays. That is mostly because I don't fully understand the broadcasting code and which code path is hit when.

[rafaqz]

Yeah. There may be another way to do this too, and we dont need it to work on more than zero dim arrays because we dont call similar for broadcasts, they stay lazy.

So the question is is this the most elegant way to handle the currently broken 0-dim case, or maybe does that matter.

Probably limiting these methods to the 0-dim edge case is safest?

[felixcremer]

We could also overload broadcast_preserve_zero_d and remove the last line of the fallback implementation because the result of Broadcast.materialize is already wrapped which it isn't in the Base Array case.

This would ensure, that a - a is always a BroadcastDiskArray also for the 0-dim case.

[sethaxen]

I'm reluctant to overload broadcast_preserve_zero_d as it's an internal function. I think just changing the 0-dim edge case would be a fine compromise.

[meggart]

Seems mergeable to me as well although I also understand @felixcremer s concerns about possible confusions in case downstream code relies on the fact the broadcast operations always return a BroadcastDiskArray. Personally I dont see too many issues, mainly for the following reasons:

1. It looks like downstream packages have not hit this code path until now, so changing the behavior probably does not break existing working code
2. We could also argue, that it is not documented, which type of struct a broadcast operation should return and that we have some freedom for changes here. In general not all operations on DIskArrays return DIskArrays (e.g. reducedim), and In general I would even be open to more flexible behavior here, e.g. even deciding to return Arrays instead of DiskArrays based the size of the return array, not only on dimensionality, although this would add some type instability. So I would suggest users can expect that whatever is returned by an operation in DiskArrays follows the AbstractArray interface but they should not expect a certain Array type. Maybe we should document this somewhere.

@sethaxen would you be up for creating a PR with your changes or would you prefer one of us adding your fixes?

[sethaxen]

@sethaxen would you be up for creating a PR with your changes or would you prefer one of us adding your fixes?

I'm happy to!

[sethaxen]

@rafaqz @felixcremer I've opened a PR fixing this issue (#286) and it's ready for review.