ND-`conv` with mixed offset/non-offset axes

[martinholters]

We currently have

julia> conv([1 1; 1 1], [1 1; 1 1])
3×3 Matrix{Int64}:
 1  2  1
 2  4  2
 1  2  1

julia> x = OffsetArray([1 1; 1 1], (-1, -1))
2×2 OffsetArray(::Matrix{Int64}, 0:1, 0:1) with eltype Int64 with indices 0:1×0:1:
 1  1
 1  1

julia> conv(x, x)
3×3 OffsetArray(::Matrix{Int64}, 0:2, 0:2) with eltype Int64 with indices 0:2×0:2:
 1  2  1
 2  4  2
 1  2  1

julia> x = OffsetArray([1 1; 1 1])
2×2 OffsetArray(::Matrix{Int64}, 1:2, 1:2) with eltype Int64 with indices 1:2×1:2:
 1  1
 1  1

julia> conv(x, x)
3×3 OffsetArray(::Matrix{Int64}, 2:4, 2:4) with eltype Int64 with indices 2:4×2:4:
 1  2  1
 2  4  2
 1  2  1

Note that the first two cases each make perfect sense on their own but are inconsistent with each other, so that conv([1 1; 1 1], [1 1; 1 1]) != conv(OffsetArray([1 1; 1 1]), OffsetArray([1 1; 1 1])) (due to different axes) as seen by that last case. Despite this inconsistency, this was deemed the least bad behavior we could think of.

The present implementation decides the indexing (with/without offset) individually per axes. For an x such that

axes(x) == (Base.OneTo(2), OffsetArrays.IdOffsetRange(values=0:1, indices=0:1))

we would (theoretically*⁾) get

axes(conv(x,x)) == (Base.OneTo(3), OffsetArrays.IdOffsetRange(values=0:2, indices=0:2))

As far as I know, there is presently no type that supports such heterogeneous axes, so it is hard to say whether this behavior is desirable. Alternatives would be

axes(conv(x,x)) == (OffsetArrays.IdOffsetRange(values=0:2, indices=0:2), OffsetArrays.IdOffsetRange(values=0:2, indices=0:2))

or

axes(conv(x,x)) == (OffsetArrays.IdOffsetRange(values=2:4, indices=2:4), OffsetArrays.IdOffsetRange(values=0:2, indices=0:2))

Without this case even being possible, I find it very hard to decide what the outcome should be. Nevertheless, we should consider what should happen if this becomes possible in the future. I tend to throw an error for now so that we are free to decide later when (and if) such an array type becomes reality.

*) Actually, with conv using similar to generate the output array and

julia> axes(similar(zeros(2,2), (0:2, Base.OneTo(3))))
(OffsetArrays.IdOffsetRange(values=0:2, indices=0:2), OffsetArrays.IdOffsetRange(values=1:3, indices=1:3))

that would be what we actually get, unless the fictitious input array type also comes with appropriate similar methods to produce a matching heterogeneous-axes array here.

Ref. #580 (comment) and the following comments, CC @wheeheee.

[martinholters]

One argument in favor of

axes(conv(x,x)) == (OffsetArrays.IdOffsetRange(values=2:4, indices=2:4), OffsetArrays.IdOffsetRange(values=0:2, indices=0:2))

might be that currently, if one of the two arrays has offset axes, we treat the other one as if it had offset axes with offset zero, so that

julia> conv([1, 1], OffsetArray([1, 1]))
3-element OffsetArray(::Vector{Int64}, 2:4) with eltype Int64 with indices 2:4:
 1
 2
 1

(and similar for more dimensions). Admittedly, I'm not convinced this is good idea either. Is this -- especially the output offset -- what one would expect?

[wheeheee]

I'm of the opinion that it's better to treat heterogeneous axes as offset by default if any are, since as you've found similar expects that too. But maybe it's easier to error for heterogeneous axes yet leave some flexibility so that the output indices calculation can be extended? As heterogeneous axes are niche anyway, whoever uses them with conv can define a new axis type and extend the functions. Something like adding

conv_with_offset(t::NTuple) = any(conv_with_offset, t)
conv_with_offset(au, av) = conv_with_offset(au) || conv_with_offset(av)

function calc_output_inds(ao, au, av)
    input_has_offset = conv_with_offset(au, av)
    if !all(input_has_offset .== conv_with_offset(ao))
        throw(ArgumentError("output must have offset axes if and only if the input has"))
    end
    offset = .!input_has_offset
    oinds = @. UnitRange{Int}(first(au) + first(av) - offset, last(au) + last(av) - offset)
    return CartesianIndices(oinds)
end

Admittedly, I'm not convinced this is good idea either. Is this -- especially the output offset -- what one would expect?

I think if this combination of arguments is allowed by default, how it works currently is reasonable enough, because the output must have offset axes, so it makes sense to treat the normal array as an offset one. And given that the 1-based behaviour is more unnatural, because Array arguments have to return Arrays, offsetting only when both axes have to be offset keeps things slightly more consistent.
Of course, it could error, but using offset axes is already an "at your own risk" kind of thing.

[martinholters]

Hm, regarding the mixing of offset-ness among inputs, I wonder whether the current behavior is what we want. Arguably, in conv, we treat arrays with non-offset axes as if they had offset -1, i.e. were zero-based. At least in the homogeneous case. Generally, one could say: "Convolution with a non-offset-axes array preserves the start index of the other array". Doesn't sound too bad. So should the offset in the mixed case be different?

julia> conv([1, 1], OffsetArray([1, 1]))
3-element OffsetArray(::Vector{Int64}, 2:4) with eltype Int64 with indices 2:4: # current
3-element OffsetArray(::Vector{Int64}, 1:3) with eltype Int64 with indices 1:3: # possible alternative
 1
 2
 1

To me, both alternatives seem equally reasonable, which makes we wonder even more whether that should just be disallowed. At least for now until a use-case manifests itself which clearly indicates which choice is to be preferred.

Another question is whether the dispatch hierarchy should be more complex to allow for more extension points, like proposed in #583 (comment). However, that seems a bit speculative and it might introduce complexity that is never used. I deliberately did not make the current functions public to allow future changes here to be non-breaking, so we can always add the complexity later when the need arises.

[martinholters]

We still might want to decide to do something useful here, but with #586, the unclear cases no throw an error, so changes in this regard won't be breaking and can be done anytime, so taking off the milestone.

[wheeheee]

Hm, regarding the mixing of offset-ness among inputs, I wonder whether the current behavior is what we want.

I've come around to disallowing mixing among inputs (permanently), after #586. One benefit is that conv would then be associative (up to numerical error), as expected, whereas if mixing is allowed, it isn't necessarily (because of different axes).