Views with IdOffsetRange as axes don't produce expected results

[jishnub]

julia> using OffsetArrays

julia> a12 = zeros(3:8, 3:4)
6×2 OffsetArray(::Array{Float64,2}, 3:8, 3:4) with eltype Float64 with indices 3:8×3:4:
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  0.0

julia> r = OffsetArrays.IdOffsetRange(Base.OneTo(3), 5)
6:8

julia> a12[r, 4] .= 3
3-element view(OffsetArray(::Array{Float64,2}, 3:8, 3:4), 6:8, 4) with eltype Float64 with indices 6:8:
 3.0
 6.9002422153289e-310
 6.9002384575573e-310

julia> a12
6×2 OffsetArray(::Array{Float64,2}, 3:8, 3:4) with eltype Float64 with indices 3:8×3:4:
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  0.0

# Expected result
julia> a12[UnitRange(r), 4] .= 3
3-element view(OffsetArray(::Array{Float64,2}, 3:8, 3:4), 6:8, 4) with eltype Float64:
 3.0
 3.0
 3.0

julia> a12
6×2 OffsetArray(::Array{Float64,2}, 3:8, 3:4) with eltype Float64 with indices 3:8×3:4:
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  3.0
 0.0  3.0
 0.0  3.0

[jishnub]

This appears to be linked to views, so it's not a broadcasting issue fundamentally.

julia> a12
6×2 OffsetArray(::Array{Float64,2}, 3:8, 3:4) with eltype Float64 with indices 3:8×3:4:
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  3.0
 0.0  3.0
 0.0  3.0

julia> v = @view a12[r, 4]
3-element view(OffsetArray(::Array{Float64,2}, 3:8, 3:4), 6:8, 4) with eltype Float64 with indices 6:8:
 3.0
 6.90023845756244e-310
 6.9002384575573e-310

julia> v2 = @view a12[UnitRange(r), 4]
3-element view(OffsetArray(::Array{Float64,2}, 3:8, 3:4), 6:8, 4) with eltype Float64:
 3.0
 3.0
 3.0

[jishnub]

Appears to be related to Base.compute_offset1.

julia> I1 = (r,4);

julia> I2 = (UnitRange(r), 4);

julia> Base.compute_offset1(a12, 1, Base.find_extended_dims(1, I1...), Base.find_extended_inds(I1...), I1)
7

julia> Base.compute_offset1(a12, 1, Base.find_extended_dims(1, I2...), Base.find_extended_inds(I2...), I2)
9

julia> @which Base.compute_offset1(a12, 1, Base.find_extended_dims(1, I2...), Base.find_extended_inds(I2...), I2)
compute_offset1(parent, stride1::Integer, dims, inds, I::Tuple) in Base at subarray.jl:377

julia> @which Base.compute_offset1(a12, 1, Base.find_extended_dims(1, I1...), Base.find_extended_inds(I1...), I1)
compute_offset1(parent, stride1::Integer, dims::Tuple{Int64}, inds::Tuple{OffsetArrays.IdOffsetRange}, I::Tuple) in OffsetArrays at /home/jishnu/.julia/packages/OffsetArrays/46NUD/src/axes.jl:133

[jishnub]

Appears to be similar to #100 , wasn't this fixed by #101 ?

[johnnychen94]

Possibly, but I'm not very confident that we should fix this in Base.compute_offset1, which IIUC is designed to calculate the offset for linear indexing. (it is possible that #101 is not a proper fix to #100)

julia> a12[r, :]
3×2 OffsetArray(::Matrix{Float64}, 6:8, 3:4) with eltype Float64 with indices 6:8×3:4:
 0.0  0.0
 0.0  0.0
 0.0  0.0

julia> a12[r, :] .= 3
3×2 view(OffsetArray(::Matrix{Float64}, 3:8, 3:4), 6:8, :) with eltype Float64 with indices 6:8×3:4:
 3.0  3.0
 3.0  3.0
 3.0  3.0

julia> a12
6×2 OffsetArray(::Matrix{Float64}, 3:8, 3:4) with eltype Float64 with indices 3:8×3:4:
 0.0  0.0
 0.0  0.0
 0.0  0.0
 3.0  3.0
 3.0  3.0
 3.0  3.0

If we think about this from another aspect, it is because that OffsetArrays has no information about what 4 represents; it could be a linear index or an offset index.

julia> r = OffsetArrays.IdOffsetRange(Base.OneTo(3), 5)
6:8

julia> a12[r, 4:4] .= 3
3×1 view(OffsetArray(::Matrix{Float64}, 3:8, 3:4), 6:8, 4:4) with eltype Float64 with indices 6:8×Base.OneTo(1):
 3.0
 3.0
 3.0

julia> a12
6×2 OffsetArray(::Matrix{Float64}, 3:8, 3:4) with eltype Float64 with indices 3:8×3:4:
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  3.0
 0.0  3.0
 0.0  3.0

For 1d cases, OffsetArrays currently interpreted it as offset index, but for nd cases, it is interpreted as a linear index... @timholy is this an intentional design?

julia> a1 = OffsetArray(collect(1:6), 3);

julia> a1[4] # offset index
1

julia> a12 = OffsetArray(reshape(collect(1:6), 2, 3), 3, 3);

julia> a12[1] # linear index
1

[jishnub]

a12[r, :] or a12[r, 4:4] are perhaps not the right comparison, as these are not FastSubArrays that uses linear indexing.

julia> (@view a12[r,:]) isa Base.FastSubArray
false

julia> (@view a12[r,4:4]) isa Base.FastSubArray
false

julia> (@view a12[r,4]) isa Base.FastSubArray
true

SubArrays use their offset only if linear indexing is used, so the first two cases will work even if IdOffsetRanges are used as indices.

julia> IndexStyle(@view a12[r,:])
IndexCartesian()

julia> IndexStyle(@view a12[r,4:4])
IndexCartesian()

julia> IndexStyle(@view a12[r,4])
IndexLinear()

I think the fix has to be either in the offset or in the indexing, since currently this leads to out-of-bounds indexing and to occasional crashes.

julia> using OffsetArrays

julia> a12 = zeros(3:8, 3:4);

julia> r = OffsetArrays.IdOffsetRange(Base.OneTo(3), 5)
6:8

julia> a12[r, :] .= 3
3×2 view(OffsetArray(::Array{Float64,2}, 3:8, 3:4), 6:8, :) with eltype Float64 with indices 6:8×3:4:
 3.0  3.0
 3.0  3.0
 3.0  3.0

julia> a12[r, 4] .= 4

signal (11): Segmentation fault
in expression starting at none:0
jl_gc_pool_alloc at /buildworker/worker/package_linux64/build/src/gc.c:1148
jl_gc_alloc_ at /buildworker/worker/package_linux64/build/src/julia_internal.h:277 [inlined]
jl_gc_alloc at /buildworker/worker/package_linux64/build/src/gc.c:3150
_new_array_ at /buildworker/worker/package_linux64/build/src/array.c:94
jl_array_copy at /buildworker/worker/package_linux64/build/src/array.c:1153
copy at ./array.jl:371 [inlined]
[...]

[jishnub]

The indexing for FastContiguousSubArrays is defined as

function getindex(V::FastContiguousSubArray, i::Int)
    @_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds r = V.parent[V.offset1 + i]
    r
end

This implies that, at least in this case with a 2D parent, V.offset1 + i should be the corresponding linear index for the parent. For example

julia> a12 = OffsetArray(reshape(1:12, 6, 2), 3:8, 3:4);

julia> v = @view a12[UnitRange(r), 4]
3-element view(OffsetArray(reshape(::UnitRange{Int64}, 6, 2), 3:8, 3:4), 6:8, 4) with eltype Int64:
 10
 11
 12

julia> (axes(v,1) |> first) + v.offset1
10

julia> Base.compute_linindex(a12, (UnitRange(r),4))
10

However this does not seem to be the case for IdOffsetRanges

julia> v2 = @view a12[r, 4]
3-element view(OffsetArray(reshape(::UnitRange{Int64}, 6, 2), 3:8, 3:4), 6:8, 4) with eltype Int64 with indices 6:8:
 13
 14
 15

julia> (axes(v2,1) |> first) + v2.offset1
13

julia> Base.compute_linindex(a12, (r,4))
10

I would have expected (axes(v2,1) |> first) + v2.offset1 to be 10 as well. Looking into the terms individually,

julia> v2.offset1
7

julia> axes(v2,1)
6:8

The axis seems fine, so it must be the offset?

[jishnub]

I think I have a solution for this, will submit a PR soon

[jishnub]

This issue exists with Base.IdentityUnitRanges as well.

julia> using OffsetArrays

julia> a = zeros(3:8, 3:4);

julia> r = Base.IdentityUnitRange(6:8);

julia> a[r, 4] .= 3

signal (11): Segmentation fault

I suppose the fix for this should go in Base.