Is OffsetArrays.jl a poison pill?

[ChrisRackauckas]

@ToucheSir

That’s totally understandable. My ask (with ismutable/ismutable wrong for FillArrays · Issue #77 · JuliaArrays/ArrayInterface.jl · GitHub as the motivating example) is that interfaces favour being conservative and maybe asking users to install a subpackage instead of generating false positives about what they support (with or without subpackages). There will always be new array packages that haven’t yet opted into the interface, and this approach would allow us as interface consumers to avoid pirating both ArrayInterface and the array packages whenever the former falls back to a too-optimistic code path and blows up.

I agree. In the end, of this needs to go to Base, and more of the definitions need errors instead of fallbacks (one of the reasons that IMO it’s not quite ready for Base).

ismutable is a surprisingly hard case though . Most people don’t use it correctly when arrays are involved.

[jishnub]

@Paul_Soderlind

For Julia 2.0 or whatever, I would like to see something like forNEW t = 4:size(x,1); i = 2:size(x,2);x[t,i] = x[t-3,i-1];... work whether x is a traditional matrix or from OffsetArrays.

To me, this approach feels like conflating cardinal vs ordinal indexing, that is, using the value of an index vs. the ordinal rank of the index in the axis. When someone writes a[2] they may mean one of two things: (a) I want the value corresponding to the key 2 (b) I want the second value in the array, irrespective of the key. The second of these requires one to use a[begin - 1 + 2]. I have written a small package that provides a convenient syntax for this:

julia> using OrdinalIndexing, OffsetArrays

julia> a = 1:10
1:10

julia> a[3rd]
3

julia> b = OffsetArray(a, -10)
1:10 with indices -9:0

julia> b[3rd]
3

Using this, we may write the loop in the OP as

julia> x = ones(-10:-1, -10:-1);

julia> for t in axes(x, 1)[4th:end], i in axes(x, 2)[2nd:end]
           x[t, i] = 2x[t-3, i-1]
       end

julia> x
10×10 OffsetArray(::Matrix{Float64}, -10:-1, -10:-1) with eltype Float64 with indices -10:-1×-10:-1:
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0
 1.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0
 1.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0
 1.0  2.0  4.0  4.0  4.0  4.0  4.0  4.0  4.0  4.0
 1.0  2.0  4.0  4.0  4.0  4.0  4.0  4.0  4.0  4.0
 1.0  2.0  4.0  4.0  4.0  4.0  4.0  4.0  4.0  4.0
 1.0  2.0  4.0  8.0  8.0  8.0  8.0  8.0  8.0  8.0

or alternately,

julia> for t in (4:size(x,1))th, i in (2:size(x,2))th
           x[t, i] = 2x[t-3, i-1]
       end

which would lead to the same result, and this should work for both standard 1-indexed arrays and OffsetArrays.

The package is not fully composable or feature-complete as of this moment, but vanilla indexing should work.

[ParadaCarleton]

Can I just say, this package looks amazing, and could prevent a lot of bugs. It would solve a lot of other problems too, e.g. the ambiguities caused when an AxisArray or similar labelled array has integer indices, and could give a very pleasing syntax for objects like ordered dictionaries or tables with named rows.

[ParadaCarleton]

In addition to the solutions offered above, I think an extra safeguard could be to use a macro instead of a new type with completely different indices. A lot of the time, you only really want the OffsetArray syntax when working locally. But when calling other people’s code, that code might assume the usual indexing, so it’s better to avoid passing the OffsetArray to them. So the macro might work like this:

lags = 4
@offset_index -lags my_3d_tensor  # indexing operations can start from -4
max_time = lastindex(my_3d_tensor, 1)
x = Vector{Matrix{Float64}}(undef, max_time)
for time in 1:max_time
    x[time] = prod(i->my_3d_tensor[i, :, :], (-lags):0)
end

Then, i could be replaced with i+lags+1 in all indexing expressions involving my_3d_tensor, making it seem like the indices start at -lags for me. But any function calls (e.g. the matrix multiplication calls from prod) would just see a perfectly normal array.