Why do views of offsetarrays have one-based indexing if bounds are specified?

[jishnub]

The index seems to be one-based or otherwise depending on whether I explicitly specify the lower and upper bounds or leave it implicit. Why is this the case?

julia> a=zeros(-1:1);

julia> b=view(a,-1:1); axes(b)
(Base.OneTo(3),)

julia> b=view(a,:); axes(b)
(Base.IdentityUnitRange(-1:1),)

julia> @which view(a,:)
view(A::AbstractArray, I::Vararg{Any,N}) where N in Base at subarray.jl:153

julia> @which view(a,-1:1)
view(A::AbstractArray, I::Vararg{Any,N}) where N in Base at subarray.jl:153

I would have expected the behavior to be identical. Why is there a difference? Is there a way to retain the index information in the view while specifying the bounds?

[jishnub]

Evidently this is how OffsetArrays are sliced, as a similar result is obtained without views.

julia> a=zeros(-1:1)
OffsetArray(::Array{Float64,1}, -1:1) with eltype Float64 with indices -1:1:
 0.0
 0.0
 0.0

julia> axes(a[:])
(Base.IdentityUnitRange(-1:1),)

julia> axes(a[-1:1])
(Base.OneTo(3),)

[tim.holy]

Indexing in Julia follows† the following rule: if b = a[idx] for some vector index idx, then b[j] = a[idx[j]]. Since -1:1 has axes of Base.OneTo(3), b has the same indexes. This is How It Should Be.

But you can do the following:

julia> using OffsetArrays

julia> a=zeros(-1:1);

julia> b = a[Base.IdentityUnitRange(-1:1)];

julia> axes(b)
(Base.IdentityUnitRange(-1:1),)

julia> axes(a)
(Base.IdentityUnitRange(-1:1),)

---

†or should follow, any violations are a bug

[ryofurue]

@tim.holy

This is How It Should Be.

What was the reasoning behind this decision? Is there a thread to read about it in this forum?

A slice of an OffsetArray loses the offsets, is quite strange (not intuitive). Those who know the b[j] == a[idx[j]] rule can deduce the How It Should Be, but you need to invoke the rule because the result is not intuitive. Those who don’t know the rule have no clue.

So, I guess there needs to be the clear discussion of this behavior in the official documentation of OffsetArrays.

Also, using Base.IdentityUnitRange() is rather bulky. I wish there were a function, let’s call it slice(), to do that: bb = slice(a, -1:1) that preserves the indexing.

[lmiq]

@ryofurue

I wish there were a function, let’s call it slice(), to do that: bb = slice(a, -1:1) that preserves the indexing.

That could be an addition to OffsetArrays, I think it is as simple as:

julia> x = OffsetArray(1:10, -2)
1:10 with indices -1:8

julia> slice(x, range) = OffsetArray(x[range], range)
slice (generic function with 1 method)

julia> y = slice(x, 2:4)
4:6 with indices 2:4

julia> y
4:6 with indices 2:4

[jishnub]

This only works with ranges. I had proposed this some years back, but it wasn’t accepted at the time. I’m not sure this is needed in the package, frankly, given that the composition works well. One may define this easily for their own use case.

Add an offset-preserving view function

[tim.holy]

@ryofurue

What was the reasoning behind this decision? Is there a thread to read about it in this forum?

b = a[idx] for some vector index idx, then b[j] = a[idx[j]]

That’s how regular array indexing works, pre OffsetArrays: if I take b = a[[3, 5, 11]], then b[1] = a[3], right? The reason is that [3, 5, 11] has indices of its own, and they are 1:3. Thus the indices of b are 1:3, i.e., the indices of the index vector become the indices of the result.

So if you want to make arrays-with-offset indices compatible with the overall AbstractArray interface, they have to follow the same rule.

[ryofurue]

@lmiq

slice(x, range) = OffsetArray(x[range], range)

But, I don’t know Julia enough to be able to handle “:” and scalar index:

using OffsetArrays
a = OffsetArray{Float64, 3}(reshape(1:60,3,4,5), -1:1, 0:3, 2:6)
ax1 = axes(a,1)
ax2 = axes(a,2)
ax3 = axes(a,3)
jslice = 2:3
b = a[:, jslice, 4] # Doesn't preserve indexing for jslice.
bb = OffsetArray(a[:, jslice, 4], ax1, jslice)  # Re-introduce jslice.
@show typeof(bb)
@show typeof(b)
@show bb[-1,3] # Fine.
@show b[-1,3] # Error

slice(x, ranges...) = OffsetArray(x[ranges...], ranges...)
bbb = slice(a, :, 2:3, 4) # Error

What to do when you don’t know the dimensionality of the OffsetArray, when you have “:”, and when one of your index “range” is actually a scalar (matrix[:,1] is a 1D vector and so OffsetArray(matrix[irange,1], irange, 1) doesn’t work).

So far, my code remains as tedious as the above example.

[ryofurue]

@tim.holy

So if you want to make arrays-with-offset indices compatible with the overall AbstractArray interface, they have to follow the same rule.

Does that mean that a lot of existing code which deals with AbstractArray already relies on the b[j] == a[idx[j]] rule? in such a way that if OffsetArrays didn’t obey the rule, significant existing library code would stop working correctly with them?

@tim.holy

That’s how regular array indexing works, pre OffsetArrays: if I take b = a[[3, 5, 11]], then b[1] = a[3], right? The reason is that [3, 5, 11] has indices of its own, and they are 1:3. Thus the indices of b are 1:3, i.e., the indices of the index vector become the indices of the result.

Thanks, but I did understand that. I wasn’t asking how the rule makes sense for 1-base arrays. It does makes sense for 1-base arrays. My point is that the rule leads to the unintuitive indexing behavior for OffsetArrays.

All in all, I guess that we should always use IdentityUnitRange with OffsetArrays. We don’t do so only because we don’t have a convenient syntax to create an IdentityUnitRange.

[bertschi]

Considering vectors as functions on finite domains, i.e., mapping indices to values, this is also just function composition:

b = a ∘ idx
# implies that
b(x) = a(idx(x))

Indeed, keys and values provide the domain and co-domain of a vector:

julia> using OffsetArrays

julia> a = zeros(-1:1);

julia> showmapping(v) = keys(v) .=> values(v)
showmapping (generic function with 1 method)

julia> showmapping(a)
3-element OffsetArray(::Vector{Pair{Int64, Float64}}, -1:1) with eltype Pair{Int64, Float64} with indices -1:1:
 -1 => 0.0
  0 => 0.0
  1 => 0.0

julia> showmapping(-1:1)
3-element Vector{Pair{Int64, Int64}}:
 1 => -1
 2 => 0
 3 => 1

julia> showmapping(Base.IdentityUnitRange(-1:1))
3-element OffsetArray(::Vector{Pair{Int64, Int64}}, -1:1) with eltype Pair{Int64, Int64} with indices -1:1:
 -1 => -1
  0 => 0
  1 => 1

# and finally the composition
julia> showmapping(a[-1:1])
3-element Vector{Pair{Int64, Float64}}:
 1 => 0.0
 2 => 0.0
 3 => 0.0

[ryofurue]

@jishnub

the composition works well

Would you care to elaborate on this? I just want to do the “composition”.

@jishnub

One may define this easily for their own use case.

I’d appreciate it if you could show how to define such view function for my cases above? Specifically, how does one define an offset_view function that can handle

bb = offset_view(a, :, 2:3, 4)

This is my real situation. Since I have’t been able to write such a convenience function myself, my current code is bulky each time I need to take a slice of a multi-dimensional OffestArray.

[jishnub]

Depending on the use case, you may use something like

julia> A = OffsetArray(rand(2,2), 1,1)
2×2 OffsetArray(::Matrix{Float64}, 2:3, 2:3) with eltype Float64 with indices 2:3×2:3:
 0.0829426  0.531739
 0.110745   0.966358

julia> function offset_view(A, inds...)
               view(A, map(to_offset_range, to_indices(A, inds))...)
       end
offset_view (generic function with 1 method)

julia> to_offset_range(r::AbstractUnitRange) = Base.IdentityUnitRange(r)
to_offset_range (generic function with 1 method)

julia> to_offset_range(x) = x
to_offset_range (generic function with 2 methods)

julia> offset_view(A, 2:3, 2:3)
2×2 view(OffsetArray(::Matrix{Float64}, 2:3, 2:3), Base.IdentityUnitRange(2:3), Base.IdentityUnitRange(2:3)) with eltype Float64 with indices 2:3×2:3:
 0.0829426  0.531739
 0.110745   0.966358

julia> offset_view(A, 2:3, 3)
2-element view(OffsetArray(::Matrix{Float64}, 2:3, 2:3), Base.IdentityUnitRange(2:3), 3) with eltype Float64 with indices 2:3:
 0.5317390043544079
 0.9663578290747614

julia> offset_view(A, 2:3, :)
2×2 view(OffsetArray(::Matrix{Float64}, 2:3, 2:3), Base.IdentityUnitRange(2:3), Base.IdentityUnitRange(Base.Slice(OffsetArrays.IdOffsetRange(values=2:3, indices=2:3)))) with eltype Float64 with indices 2:3×Base.Slice(OffsetArrays.IdOffsetRange(values=2:3, indices=2:3)):
 0.0829426  0.531739
 0.110745   0.966358

[Tamas_Papp]

I don’t think OffsetArrays supports view directly, but you could just wrap them, as in

using OffsetArrays

_drop(::Integer, ι...) = (_drop(ι...)...,)
_drop(i1, ι...) = (i1, _drop(ι...)...,)

function offset_view(a::AbstractArray, indices...)
    OffsetArray(view(a, indices...), _drop(indices...)...)
end

The _drop is required to make sure integer indices are dropped. I did not test this extensively, so there may be horrible corner cases.

[ryofurue]

@jishnub

Depending on the use case, you may use something like . . .

Thanks! But your code indicates that such a function should be provided by the OffestArrays package because your code is too esoteric for an ordinary user to come up with.

So, I disagree with your reasoning that such a function isn’t a useful part of the package because the user can easily write it.

I’m sure that I will not be able to write such a function “Depending on [my] use case”.

[mbauman]

I mean, it’s just a very different meaning of a view (and indexing in general).

The “consistent” behavior would mean that A[2:3] should return an offset array with indices 2:3, regardless of the type or indices of A. That’s just a different sort of indexing behavior than what Julia has standardized on.