juliaarrays/offsetarrays.jl

Fortran-like arrays with arbitrary, zero or negative starting indices.

offsetarrays.jl's Issues

`searchsorted` causes segfault

julia> using OffsetArrays

julia> o = OffsetArray([2,2,3],-1:1)
3-element OffsetArray(::Array{Int64,1}, -1:1) with eltype Int64 with indices -1:1:
 
2
 2
 3

julia> searchsorted(o,1)
Segmentation fault: 11

Using the debugger reveals its a strange bug:

julia> @enter searchsorted(o,1)
In #searchsorted#6(lt, by, rev, order, , v, x) at sort.jl:296
>1  1 ─ %1 = (Base.Order.ord)(lt, by, rev, order)
 2  │   %2 = (searchsorted)(v, x, %1)
 3  └──      return %2

About to run: (Base.Order.ord)(isless, identity, nothing, Base.Order.ForwardOrdering())
1|debug> n
ERROR: BoundsError: attempt to access 3-element OffsetArray(::Array{Int64,1}, -1:1) with eltype Int64 with indices -1:1 at index [9223372036854775807]
Stacktrace:
 [1] throw_boundserror(::OffsetArray{Int64,1,Array{Int64,1}}, ::Tuple{Int64}) at abstractarray.jl:538
 [2] checkbounds(::OffsetArray{Int64,1,Array{Int64,1}}, ::Tuple{Int64}) at abstractarray.jl:503
 [3] getindex(::OffsetArray{Int64,1,Array{Int64,1}}, ::Int64) at /Users/solver/.julia/packages/OffsetArrays/vIbpP/src/OffsetArrays.jl:135
 [4] searchsorted(::OffsetArray{Int64,1,Array{Int64,1}}, ::Int64, ::Int64, ::Int64, ::Base.Order.ForwardOrdering) at sort.jl:214
 [5] searchsorted(::OffsetArray{Int64,1,Array{Int64,1}}, ::Int64, ::Base.Order.ForwardOrdering) at sort.jl:294
 [6] #searchsorted#6(::typeof(isless), ::typeof(identity), ::Nothing, ::Base.Order.ForwardOrdering, ::typeof(searchsorted), ::OffsetArray{Int64,1,Array{Int64,1}}, ::Int64) at sort.jl:296
 [7] searchsorted(::OffsetArray{Int64,1,Array{Int64,1}}, ::Int64) at sort.jl:296

PS I want this to support offset BlockArrays: JuliaArrays/BlockArrays.jl#95

collapse nested SubArray and OffsetArray

import OffsetArrays: OffsetArray

a = rand(4,4)
a[:] .= 1:16

sa = view(a, 1:3, 1:3)
osa = OffsetArray(sa, (11:13, 11:13))
sosa = view(osa, 11:12, 11:12)

function Base.view(osa::OffsetArray{T, N, <:SubArray}, I::Vararg{Any,N}) where {T, N}
    view(parent(osa), map((a,b) -> a .- b, I, osa.offsets)...)
end

sosa2 = view(osa, 11:12, 11:12)

julia> sosa
2×2 view(OffsetArray(view(::Array{Float64,2}, 1:3, 1:3), 11:13, 11:13), 11:12, 11:12) with eltype Float64:
 1.0  5.0
 2.0  6.0

julia> sosa2
2×2 view(::Array{Float64,2}, 1:2, 1:2) with eltype Float64:
 1.0  5.0
 2.0  6.0

The map and splat is probably not great for performance. Other than that, does this seem like a reasonable idea?

Should we allow non-integer offsets?

If you have the vector a = rand(5) and perform finite differencing a[i+1]-a[i], the difference values are best considered as being defined on the half-grid. I wonder whether it makes sense to define indices(da, 1) == 1.5:1:4.5 in such cases.

`show(::OffsetArrays.IdOffsetRange` should show its type

import OffsetArrays

julia> Base.OneTo(3)
Base.OneTo(3)

julia> Base.IdentityUnitRange(1:3)
Base.IdentityUnitRange(1:3)

julia> OffsetArrays.IdOffsetRange(1:3)
1:3

I think this might be misleading, as I mentioned in #120 (comment)

Reduction over OffsetArray views triggers typeassert error

On OffsetArrays master, Julia 1.3.0,

julia> sum(view(OffsetArray(rand(3,3,3), 0, 0, 0), :, :, 1:2), dims=(2,3))

ERROR: TypeError: in typeassert, expected Tuple{Base.IdentityUnitRange{UnitRange{Int64}},Base.IdentityUnitRange{UnitRange{Int64}},Base.OneTo{Int64}}, got Tuple{Base.IdentityUnitRange{UnitRange{Int64}},UnitRange{Int64},Base.OneTo{Int64}}
Stacktrace:
 [1] reduced_indices(::Tuple{Base.IdentityUnitRange{UnitRange{Int64}},Base.IdentityUnitRange{UnitRange{Int64}},Base.OneTo{Int64}}, ::Tuple{Int64,Int64}) at .\reducedim.jl:57
 [2] reduced_indices at .\reducedim.jl:15 [inlined]
 [3] reducedim_initarray(::SubArray{Float64,3,OffsetArray{Float64,3,Array{Float64,3}},Tuple{Base.Slice{Base.IdentityUnitRange{UnitRange{Int64}}},Base.Slice{Base.IdentityUnitRange{UnitRange{Int64}}},UnitRange{Int64}},true}, ::Tuple{Int64,Int64}, ::Float64, ::Type{Float64}) at .\reducedim.jl:92
 [4] reducedim_initarray at .\reducedim.jl:93 [inlined]
 [5] reducedim_init at .\reducedim.jl:172 [inlined]
 [6] _mapreduce_dim at .\reducedim.jl:317 [inlined]
 [7] #mapreduce#584 at .\reducedim.jl:307 [inlined]
 [8] #mapreduce at .\none:0 [inlined]
 [9] _sum at .\reducedim.jl:679 [inlined]

matrix multiplication not supported?

julia> A = OffsetArray([1 2; 3 4], 0:1, 0:1)
2×2 OffsetArray(::Array{Int64,2}, 0:1, 0:1) with eltype Int64 with indices 0:1×0:1:
 1  2
 3  4

julia> b = OffsetArray([5; 6], 0:1)
2-element OffsetArray(::Array{Int64,1}, 0:1) with eltype Int64 with indices 0:1:
 5
 6

julia> A*b
ERROR: ArgumentError: offset arrays are not supported but got an array with index other than 1
Stacktrace:
 [1] require_one_based_indexing at .\abstractarray.jl:89 [inlined]
 [2] generic_matvecmul!(::OffsetArray{Int64,1,Array{Int64,1}}, ::Char, ::OffsetArray{Int64,2,Array{Int64,2}}, ::OffsetArray{Int64,1,Array{Int64,1}}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\LinearAlgebra\src\matmul.jl:501
 [3] mul! at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\LinearAlgebra\src\matmul.jl:77 [inlined]
 [4] *(::OffsetArray{Int64,2,Array{Int64,2}}, ::OffsetArray{Int64,1,Array{Int64,1}}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\LinearAlgebra\src\matmul.jl:51
 [5] top-level scope at none:0

Equality of empty OffsetArrays changes in Julia 1.5

Previously, empty arrays with different offsets compared as different because their axes compared as different.

In Julia master (JuliaLang/julia#32348), empty ranges are now equal, making OffsetArray([],-1) == OffsetArray([], -2).

I'm not sure what the right semantic is here; I believe this was untested in the entire ecosystem.

Why is Diagonal(::OffsetArray) not supported?

At the moment Diagonal(::OffsetArray) is explicitly forbidden:

julia> using OffsetArrays

julia> x = OffsetArray(randn(5),-1)
5-element OffsetArray(::Array{Float64,1}, 0:4) with eltype Float64 with indices 0:4:
 -1.1571843743724437
  1.182586195865172
 -1.7939602910206534
 -0.42805091122795136
 -1.2399140299181175

julia> Diagonal(x)
ERROR: ArgumentError: offset arrays are not supported but got an array with index other than 1
Stacktrace:
 [1] require_one_based_indexing at ./abstractarray.jl:89 [inlined]
 [2] Diagonal at /Users/solver/Projects/julia-1.4/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/diagonal.jl:9 [inlined]
 [3] Diagonal(::OffsetArray{Float64,1,Array{Float64,1}}) at /Users/solver/Projects/julia-1.4/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/diagonal.jl:13
 [4] top-level scope at REPL[3]:1

However, supporting it seems to require just 2 changes in Base:

julia> Base.require_one_based_indexing(::OffsetArray) = nothing
julia> Base.axes(D::Diagonal) = (axes(D.diag,1),axes(D.diag,1))

julia> Diagonal(x)
5×5 Diagonal{Float64,OffsetArray{Float64,1,Array{Float64,1}}} with indices 0:4×0:4:
 -1.15718   ⋅         ⋅         ⋅          ⋅ 
   ⋅       1.18259    ⋅         ⋅          ⋅ 
   ⋅        ⋅       -1.79396    ⋅          ⋅ 
   ⋅        ⋅         ⋅       -0.428051    ⋅ 
   ⋅        ⋅         ⋅         ⋅        -1.23991

julia> Diagonal(x)[0,0]
-1.1571843743724437

Is there any reason this was not done?

One reason I ask is that the definition Base.axes(D::Diagonal) = (axes(D.diag,1),axes(D.diag,1)) is helpful for BlockArrays.jl, so that Diagonal(::BlockArray) automatically inherits the block structure of the diagonal. I was thinking of making a PR into Base to propose this change but wanted to know context why it wasn't done already.

SubArrays for OffsetVectors broken

Example

julia> g = OffsetArray(full(-2:3), (-2:3,))
OffsetArrays.OffsetArray{Int64,1,UnitRange{Int64}} with indices -2:3:
 -2
 -1
  0
  1
  2
  3
julia> view(g, -1:2)
4-element SubArray{Int64,1,OffsetArrays.OffsetArray{Int64,1,UnitRange{Int64}},Tuple{UnitRange{Int64}},true}:
   2
   3
 #undef
 #undef

It seems to be calculating the offset1 wrong for fast linear indexing. With quick look at the Julia code, I found a suspicious comment:

compute_offset1(parent, stride1::Integer, dims, inds, I::Tuple) =
    (@_inline_meta; compute_linindex(parent, I) - stride1)  # linear indexing starts with 1

Missing `diff` for OffsetVector

I noted Base.diff does not work for OffsetArrays. MWE

julia> x = OffsetVector([0, 1, 2, 3], 0:3)
julia> diff(x)
ERROR: ArgumentError: offset axes unsupported
Stacktrace:
 [1] #diff#395(::Int64, ::Function, ::OffsetArray{Int64,1,Array{Int64,1}}) at ./multidimensional.jl:741
 [2] #diff at ./none:0 [inlined]
 [3] diff(::OffsetArray{Int64,1,Array{Int64,1}}) at ./multidimensional.jl:708
 [4] top-level scope at none:0

It would be nice to have this.

Allow for constructing arrays with CartesianRange

Something like

OffsetArray{Float64}(undef, CartesianIndex(3,4,-1):CartesianIndex(7,9,5))

Taking views of OffsetArrays

My package is making use of the following call:

using Images

D = imfilter(ones(10,10), centered(ones(3,3)))
view(D,1:2,:) # ERROR: indices go from 2 to 9

I am not very familiar with the concept of offset arrays, but I am using views everywhere in my code to save memory allocations: https://github.com/juliohm/ImageQuilting.jl/blob/master/src/iqsim.jl#L258-L259

In the package, the D can either be a normal Array if I use "symmetric" padding or an OffsetArray if I use Inner() padding, both defined in ImageFiltering.jl. Can you help solve this case-by-case issue, or how to convert OffsetArrays to normal Arrays?

Performance hit in computing axes of views of OffsetArrays with slices as axes

julia> a = ones(3,5);

julia> oa = OffsetArray(a, axes(a));

julia> voa = @view oa[:,:];

julia> @btime axes($(Ref(voa))[]);
  14.806 ns (0 allocations: 0 bytes

Compared to this:

julia> voa = @view oa[1:2,1:2];

julia> @btime axes($(Ref(voa))[]);
  3.740 ns (0 allocations: 0 bytes)

julia> voa = @view oa[:,1:2];

julia> @btime axes($(Ref(voa))[]);
  3.213 ns (0 allocations: 0 bytes)

julia> voa = @view oa[1:2,:];

julia> @btime axes($(Ref(voa))[]);
  3.211 ns (0 allocations: 0 bytes)

Why is the axes computation so much more expensive when both the axes are slices? Note that this is not the case with Arrays:

julia> va = @view a[:,:];

julia> @btime axes($(Ref(va))[]);
  3.203 ns (0 allocations: 0 bytes)

Concatenation of 3D OffsetArrays fails if their axes do not start at 1

I try to hcat an OffsetArray with indices (0:-1, 0:-1,0:1) and size=(0,0,2) with another OffsetArray with indices (0:-1,0:15,0:1) and size(0,16,2) and I get the error ERROR: BoundsError: attempt to access 0×16×2 Array{Int8,3} at index [1:0, 0:15, 1:2]

julia> form=OffsetArray(reshape(zeros(Int8,0),0,0,2),0:-1,0:-1,0:1)
julia> exp=OffsetArray(reshape(zeros(Int8,0),0,16,2),0:-1,0:15,0:1)
julia> hcat(form,exp) Error!!

julia> hcat(form,exp)
ERROR: BoundsError: attempt to access 0×16×2 Array{Int8,3} at index [1:0, 0:15, 1:2]
Stacktrace:
 [1] throw_boundserror(::Array{Int8,3}, ::Tuple{UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}}) at ./abstractarray.jl:484
 [2] checkbounds at ./abstractarray.jl:449 [inlined]
 [3] _setindex! at ./multidimensional.jl:638 [inlined]
 [4] setindex! at ./abstractarray.jl:998 [inlined]
 [5] __cat(::Array{Int8,3}, ::Tuple{Int64,Int64,Int64}, ::Tuple{Bool,Bool}, ::OffsetArray{Int8,3,Array{Int8,3}}, ::Vararg{OffsetArray{Int8,3,Array{Int8,3}},N} where N) at ./abstractarray.jl:1379
 [6] _cat_t(::Val{2}, ::Type, ::OffsetArray{Int8,3,Array{Int8,3}}, ::Vararg{OffsetArray{Int8,3,Array{Int8,3}},N} where N) at ./abstractarray.jl:1361
 [7] #cat_t#99(::Val{2}, ::Function, ::Type{Int8}, ::OffsetArray{Int8,3,Array{Int8,3}}, ::Vararg{OffsetArray{Int8,3,Array{Int8,3}},N} where N) at ./abstractarray.jl:1353
 [8] (::getfield(Base, Symbol("#kw##cat_t")))(::NamedTuple{(:dims,),Tuple{Val{2}}}, ::typeof(Base.cat_t), ::Type{Int8}, ::OffsetArray{Int8,3,Array{Int8,3}}, ::Vararg{OffsetArray{Int8,3,Array{Int8,3}},N} where N) at ./none:0
 [9] _cat at ./abstractarray.jl:1481 [inlined]
 [10] #cat#100 at ./abstractarray.jl:1480 [inlined]
 [11] #cat at ./none:0 [inlined]
 [12] hcat(::OffsetArray{Int8,3,Array{Int8,3}}, ::OffsetArray{Int8,3,Array{Int8,3}}) at ./abstractarray.jl:1489
 [13] top-level scope at none:0

However if I try with arrays with regular indices

julia> form=reshape(zeros(Int8,0),0,0,2)
julia> exp=reshape(zeros(Int8,0),0,16,2)
juilia> hcat(form,exp) works!!!

Thanks

Reshape with Colon doesn't work, can segfault

The reshape call with a Colon doesn't seem to be collected by OffsetArrays, thus due to another bug in Base, one can get segfaults.
See /JuliaLang/julia/issues/33603 and /JuliaLang/julia/issues/33614

bad printing style for column matrices in tuples (probably unintended)

As you can see we have some duplication of the semicolons:

julia> (ones(3,1),)
([1.0; 1.0; 1.0],)

julia> (centered(ones(3,1)),)
([1.0; ; 1.0; ; 1.0],)

EDIT: doing println also doesn't work

julia> println(centered(ones(3,1)))
[1.0; ; 1.0; ; 1.0]

julia> println(ones(3,1))
[1.0; 1.0; 1.0]

Offset map

I have a very simple suggestion, based on how I almost always use OffsetArrays. Would there be any interest adding the utility function:

offsetmap(f,ind) = OffsetArray(map(f,ind),ind)

support offsett-version of linear index?

I suppose we need to support offset-version of linear indexing described in https://docs.julialang.org/en/v1/manual/arrays/#Omitted-and-extra-indices-1 ?

julia> m = OffsetArray(rand(3, 3, 1), -2, -2, -1)
OffsetArray(::Array{Float64,3}, -1:1, -1:1, 0:0) with eltype Float64 with indices -1:1×-1:1×0:0:
[:, :, 0] =
 0.087044  0.377619  0.535544
 0.363213  0.84182   0.440447
 0.134247  0.849669  0.561989

julia> m[0,0] # omit trailing length-1 dimensions 😄 
0.8418204286682247

julia> m[0,0,0]
0.8418204286682247

julia> m[0] # not-offsetted linear index 😢 
ERROR: BoundsError: attempt to access OffsetArray(::Array{Float64,3}, -1:1, -1:1, 0:0) with eltype Float64 with indices -1:1×-1:1×0:0 at index [0]
Stacktrace:
 [1] throw_boundserror(::OffsetArray{Float64,3,Array{Float64,3}}, ::Tuple{Int64}) at ./abstractarray.jl:484
 [2] checkbounds at ./abstractarray.jl:449 [inlined]
 [3] getindex(::OffsetArray{Float64,3,Array{Float64,3}}, ::Int64) at /Users/jc/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:140
 [4] top-level scope at none:0

julia> m[5]
0.8418204286682247

Constrain the types dispatched on in Base functions like fill, ones, zeros, trues, falses

Currently the function signature for fill is

OffsetArrays.jl/src/OffsetArrays.jl

Lines 212 to 213 in 7294f3d

	Base.fill(x, inds::Tuple{UnitRange,Vararg{UnitRange}}) =
	fill!(OffsetArray{typeof(x)}(undef, inds), x)

However the function clearly only works with integral values, as the OffsetArray constructor only accepts Ints. Something like this will fail:

julia> fill(3.0,UnitRange(3.0,5.0))
ERROR: MethodError: no method matching OffsetArray(::Array{Float64,1}, ::Tuple{Float64})

Should the types dispatched upon be constrained to UnitRange{<:Integer}? This will also allow other packages to define these methods for various other UnitRange types.

The same for other functions like zeros, ones, trues, falses.

Performance of OffsetArrays

Hello,

I have posted on discourse but it seems to have gone unnoticed. Apologies if it is inappropriate to file multiple posts. Cross-referencing here:
https://discourse.julialang.org/t/performance-of-offsetarrays/5769?u=raminammour

Clarify `==`and `===` for `axes` tests

In the tests I see that there are some places where === is used to compare results of axes, which is good since 1:3 == Base.OneTo(3), and some places want to test specifically for Base.OneTo(3)

But I also see many cases of e.g. a == IdOffsetRange(3:5) and a == IdentityUnitRange(-4:-3), and I don't see the point of doing that instead of a == 3:5, unless the test is really intended to check for ===, in which case in at least a few places it should be a === OffsetArrays.IdOffsetRange(Base.OneTo(3), 2).

For example:

OffsetArrays.jl/test/runtests.jl

Line 446 in 7294f3d

@test axes(v) == (Base.OneTo(1), IdentityUnitRange(-4:-3))

Some questions on useability

Well, firstly thanks for providing OffsetArrays.

I want to make it simpler/cleaner to specify OffsetArrays backed by standard arrays. And also use analogous syntax to standard Arrays for creating and typing OffsetArrays.

An example might be best:

type MyType
    x::OffsetStdVector{Int} 
    y::OffsetStdArray{Float64,2}
end

a = OffsetVector{Int}(0:5) 
b = OffsetArray{Float64}(0:12, -1:1) 
c = MyType(a, b)

a = zeros(Int, 0:5)    # [3]
b = zeros(Float64, 0:12, -1:1)    #[3]
c = MyType(a, b)

This can be achieved currently with something like:

using OffsetArrays
import OffsetArrays.OffsetVector

typealias OffsetStdVector{T}  OffsetArray{T,1,Array{T,1}}
typealias OffsetStdArray{T,N} OffsetArray{T,N,Array{T,N}}

(::Type{OffsetVector{T}}){T}(inds::UnitRange{Int}) = OffsetArray{T}(inds)
(::Type{OffsetArray{T}}){T}(inds::UnitRange{Int}...) = OffsetArray{T}(inds)

Base.zeros{N}(T::Type, inds::Base.Indices{N}) = fill!(OffsetArray{T}(inds), zero(T))
Base.zeros(T::Type, inds::UnitRange{Int}...) = zeros(T, inds)

Ideally OffsetArray could be renamed to say OffsetGenArray, so that standard array-backed type could be just OffsetArray

In any case,

1. Could OffsetVector be exported?
2. Is there any problem with defining methods so that OffsetArray{T}(inds) can construct (similar to recommended way to construct standard arrays)?
3. Why is zero deprecated? using syntax as above seems quite natural to me.

Printing with :limit = true is broken

julia> show(IOContext(stdout, :limit=>true), OffsetArray(collect(0:20), 0:20))
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  11, 12, 13, 14, 15, 16, 17, 18, 19, 20]
julia> OffsetArray(collect(0:20), 0:20)
OffsetArray(::Array{Int64,1}, 0:20) with eltype Int64 with indices 0:20:
  0
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20

(note that the first entry in the limited print is 1, not 0 as expected).

Broadcasting with IdOffsetRange seems wrong

This seems fine, although the error message is a bit odd:

oa = OffsetArray(1:4,1)
axes(oa .+ 1)  # (2:5)
oa .+ (1:4)   # 
DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 4 and 4")

But:

oa_axes = axes(oa)[1]  # 2:5
axes(oa_axes)  # (2:5)
axes(oa_axes .+ 1)  # (Base.OneTo(4))
oa_axes .+ (1:4)  #  3:2:9

I would expect oa_axes to behaves like oa in broadcasting.

support CartesianIndex offset

The following usage could be helpful to deal with multi-dimensional arrays

OffsetArray(rand(4, 4), CartesianIndex(-1, -1))

I could do it later if you think this is a good idea.

Iterating using `axes` causes a lot of memory allocations

I upgraded from 0.11.4 to 1.0.1 and noticed that my code slowed down significantly. Here is a MWE

using Pkg
using OffsetArrays

println(Pkg.installed()["OffsetArrays"])

function foo(X)
    s = 0.
    for j in 1:1024, i in 1:1024
        for xj in axes(X, 2), xi in axes(X, 1)
            s += X[xi, xj]
        end
    end
    s
end

X = OffsetArray(randn(5, 5), -2:2, -2:2)
foo(X)
@time foo(X)

Results:

0.11.4
  0.036787 seconds (5 allocations: 176 bytes)

1.0.0
  1.207429 seconds (56.62 M allocations: 1.688 GiB, 13.36% gc time)

1.0.1
  1.230504 seconds (56.62 M allocations: 1.688 GiB, 13.13% gc time)

I can resort to firstindex, lastindex for now.

`setindex_shape_check` error

using OffsetArrays

a = OffsetArray{Complex{Float64}}(undef, 100, -5:5, 1:3)
a[:, 0, :] # works
a[:, 0, :] = 3.0

gives

ERROR: MethodError: no method matching setindex_shape_check(::Float64, ::Int64, ::Int64, ::Int64)
Closest candidates are:
  setindex_shape_check(::AbstractArray, ::Integer...) at indices.jl:154
  setindex_shape_check(::AbstractArray{#s57,1} where #s57, ::Integer, ::Integer) at indices.jl:196
  setindex_shape_check(::AbstractArray{#s57,2} where #s57, ::Integer, ::Integer) at indices.jl:200
  ...
Stacktrace:
 [1] macro expansion at ./multidimensional.jl:641 [inlined]
 [2] _unsafe_setindex!(::IndexLinear, ::OffsetArray{Complex{Float64},3,Array{Complex{Float64},3}}, ::Float64, ::Base.Slice{UnitRange{Int64}}, ::Int64, ::Base.Slice{UnitRange{Int64}}) at ./multidimensional.jl:636
 [3] _setindex! at ./multidimensional.jl:631 [inlined]
 [4] setindex!(::OffsetArray{Complex{Float64},3,Array{Complex{Float64},3}}, ::Float64, ::Function, ::Int64, ::Function) at ./abstractarray.jl:998
 [5] top-level scope at none:0

How can I avoid this error? Am I doing something I shouldn't?

Julia 1.0, Ubuntu 16.04

BTW thanks immensely for this package; I'm writing Ewald sums and these would be difficult without the -5:5 array queries.

LinearIndexing for a zero-dimensional view of an OffsetVector doesn't work as expected

For a zero-dimensional array a we usually expect a[] == a[1]. However this does not work with zero-D views of OffsetVectors

julia> a = OffsetArray(1:3, 0:2);

julia> b = @view a[0]
0-dimensional view(OffsetArray(::UnitRange{Int64}, 0:2), 0) with eltype Int64:
1

julia> b[]
1

julia> b[1]
2

This happens because for FastContiguousSubArrays, the indexing produces b[1] -> a[1]. Perhaps the offset calculation for the view needs to be different for OffsetVectors?

It appears adding the method

Base.compute_offset1(parent::OffsetVector, stride1::Integer, I::Tuple{Int}) = I[1] - stride1

fixes the issue in this particular case, however I'm not sure if this holds in general. I wonder if this is the correct resolution? In that case I can consider a PR.

Keep indices when slicing

It is often useful to take a subset of an OffsetArray like this:

using OffsetArrays
full = OffsetArray(ones(101, 101), (-50:50, -50:50))
part = full[-10:10, -10:10]
# would like to pass part to a function or somewhere else,
# and access it as e.g. part[-5, 10]

However, indexing into an OffsetArray returns a regular Array with one-based indices. To keep the offset I currently do part = OffsetArray(full[-10:10, -10:10], (-10:10, -10:10)).

I understand that the current behavior is for consistency with slicing regular Arrays. But is there a more convenient, automatic way to slice OffsetArray and keep indices properly offset? If not - I looked into code and didn't find anything - how do you think this would be implemented best?

StackOverflowError in similar with Colons

julia> oa = OffsetArray(rand(3,3),-1:1,-1:1)
3×3 OffsetArray(::Array{Float64,2}, -1:1, -1:1) with eltype Float64 with indices -1:1×-1:1:
 0.780486  0.736212  0.62076
 0.465256  0.477541  0.478559
 0.839068  0.50454   0.0118495

julia> similar(oa, Float64, (:,:))
ERROR: StackOverflowError:
Stacktrace:
 [1] similar(::OffsetArray{Float64,2,Array{Float64,2}}, ::Type{Float64}, ::Tuple{Colon,Colon}) at /home/jishnu/.julia/packages/OffsetArrays/Z45he/src/OffsetArrays.jl:102 (repeats 79984 times)

Presumably this happens because

OffsetArrays.jl/src/OffsetArrays.jl

Lines 101 to 104 in bd7608e

	function Base.similar(A::AbstractArray, ::Type{T}, inds::Tuple{OffsetAxis,Vararg{OffsetAxis}}) where T
	B = similar(A, T, map(indexlength, inds))
	return OffsetArray(B, map(offset, axes(B), inds))
	end

and

OffsetArrays.jl/src/OffsetArrays.jl

Line 225 in bd7608e

indexlength(i::Colon) = Colon()

lbound and ubound functions

Fortran arrays has lbound and ubound. Is it exist with OffsetArrays ?

I code it myself like this

lbound( array :: OffsetArray, dim ) = axes(array)[dim].indices[1]
ubound( array :: OffsetArray, dim ) = axes(array)[dim].indices[end] 

Is there a better way to do this ?

thank you very much for bringing this feature to Julia. Very practical especially for block parallelization and manage ghost points.

Wrong results using argmin or argmax

Using argmin or argmax along a dimension of an OffsetArray gives wrong results for some offsets
Examples:

using OffsetArrays
A = [1 2; 1 2]
for offset in 0:-1:-3
    A_off = OffsetArray(A, offset, -1)
    A_off[argmin(A_off, dims=2)] |> println
end

produces

[1; ; 1] # Correct
[2; ; 1] # Wrong
[1; ; 2] # Wrong
[1; ; 1] # Correct

Inference with other custom-indexing types?

I have implemented custom circular buffer type with cumulative indexing, so that it has offset field that grows each time an element is pushed back to the full buffer. (Why I did not used OffsetArray because of it's immutable and I can't change offsets tuple).

Then, I figured that I miss some similar methods, https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array-1

And that they are used from this package, automatically converting my buffer to OffsetArray on some operations, due to this method:

Base.similar(::Type{T}, shape::Tuple{OffsetAxis,Vararg{OffsetAxis}}) where {T<:AbstractArray} =
    OffsetArray(T(undef, map(indexlength, shape)), map(indexoffset, shape))

with OffsetAxis broadly defined as:

const OffsetAxis = Union{Integer, UnitRange, Base.OneTo, IdentityUnitRange}

How should I write those similar and axes methods for my buffer type to not have any inference between packages?

Views with IdOffsetRange as axes don't produce expected results

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

BoundsError on type-unstable map

Each of these lines

map(!, OffsetArray([true,missing],2))
map(x -> x>1 ? 1 : nothing, OffsetArray(1:2,2))

errors with

BoundsError: attempt to access 2-element OffsetArray(::Array{Union{Missing, Bool},1}, 3:4) with eltype Union{Missing, Bool} with indices 3:4 at index [3:5]

`diagm` doesn't work

A = OffsetArray(rand(3), -3)
diagm(A)

results in

ERROR: size not supported for arrays with indices (-2:0,); see http://docs.julialang.org/en/latest/devdocs/offset-arrays/
Stacktrace:
 [1] errmsg(::OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}) at /Users/twan/code/julia/MixedIntegerExperiments/v0.6/OffsetArrays/src/OffsetArrays.jl:48
 [2] diagm at ./linalg/dense.jl:250 [inlined] (repeats 2 times)

on nightly. This should probably be fixed in Base, but I wanted to get your input first.

As for the desired output type of diagm, I guess it should always be

typeof(similar(A, first(indices(A)), first(indices(A))))

?

OffsetArrays.jl package does not work in 0.3.12 and 0.4.1

'''
Version 0.3.12 , using console
julia> Pkg.installed()
Dict{ASCIIString,VersionNumber} with 1 entry:
"OffsetArrays" => v"0.0.0-"
julia> using OffsetArrays
julia> y = OffsetArray(Float32,-1:1);
julia> y
NaN
NaN
0.0
julia> y[0] = 5.0
ERROR: arrayset not defined
in setindex! at /home/juser/.julia/v0.3/OffsetArrays/src/OffsetArrays.jl:72
julia>

Version 0.3.12, using Ijulia
In: Pkg.installed()
Dict{ASCIIString,VersionNumber} with 1 entry:
"OffsetArrays" => v"0.0.0-"
In: using OffsetArrays
In: y = OffsetArray(Float32, -1:1)
arraysize not defined

in showarray at show.jl:1097
in anonymous at replutil.jl:27
in with_output_limit at ./show.jl:1136
in writemime at replutil.jl:26
in writemime at multimedia.jl:41
in sprint at iostream.jl:229
in display_dict at /opt/julia_packages/.julia/v0.3/IJulia/src/execute_request.jl:26
In:

Version 0.4.1, using ijulia
In: Pkg.installed()
Dict{ASCIIString,VersionNumber} with 1 entry:
"OffsetArrays" => v"0.0.0-"
In: using OffsetArrays
LoadError: LoadError: UndefVarError: Range1 not defined
while loading /home/juser/.julia/v0.4/OffsetArrays/src/OffsetArrays.jl, in expression starting on line 8
while loading In[6], in expression starting on line 1

in include at ./boot.jl:261
in include_from_node1 at ./loading.jl:304
in require at ./loading.jl:243
In:
'''

Type-piracy in similar with `OffsetAxis`

This is type-piracy:

OffsetArrays.jl/src/OffsetArrays.jl

Line 99 in d50463f

function Base.similar(A::AbstractArray, ::Type{T}, inds::Tuple{OffsetAxis,Vararg{OffsetAxis}}) where T

enumerate returns wrong indices

I tried using OffsetArrays with enumerate, but the indices are wrong. Here is a small example

x = collect(1:4)
collect(enumerate(OffsetArray(x[2:4], 2:4)))

which gives

3-element Array{Tuple{Int64,Int64},1}:
 (1, 2)
 (2, 3)
 (3, 4)

I am not sure if the problem is with enumerate or with OffsetArrays.

setindex! should return the array, not the value

Eg.

julia> setindex!(ones(3),2,1)
3-element Array{Float64,1}:
 2.0
 1.0
 1.0

julia> setindex!(OffsetArray(ones(3),1:3),2,1)
2

It might be better if setindex! returned the array to be consistent with Base.

How to add arrays with different offsets?

I wanted to add two arrays with different offsets assuming zero-padding, and tried to use paddedviews(). But, didn't work :-(
So, I tried to create a function to do that... Below is my first attempt in a one-dimentional case.
Before extending it to multi-dimensional arrays, it looks already pretty complicated :-(

Are there better ways to do that?

using OffsetArrays
a = [1,2,3]; b = [4,5,6];
a_o = OffsetArray(a, -2); b_o = OffsetArray(b, 1);
function add_offsets(a_o, b_o)
    a_r = indices(a_o); b_r = indices(b_o);
    a_min = minimum.(a_r)[1]; a_max = maximum.(a_r)[1];
    b_min = minimum.(b_r)[1]; b_max = maximum.(b_r)[1];
    ab_min = min(a_min, b_min); ab_max = max(a_max, b_max);
    a_p = zeros(eltype(a_o), ab_max - ab_min + 1)
    a_p = OffsetArray(a_p, ab_min:ab_max)
    a_p[a_min:a_max] = a_o[a_min:a_max]
    a_p[b_min:b_max] += b_o[b_min:b_max]
    a_p
end
add_offsets(a_o, b_o)

Multi-arg eachindex results in error

On both 0.5 and nightly:

x = rand(3)
y = OffsetArray(x, -length(x))
eachindex(x, y)

result on nightly:

ERROR: size not supported for arrays with indices (-2:0,); see http://docs.julialang.org/en/latest/devdocs/offset-arrays/
Stacktrace:
 [1] errmsg(::OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}) at /Users/twan/code/julia/MixedIntegerExperiments/v0.6/OffsetArrays/src/OffsetArrays.jl:48
 [2] eachindex(::Array{Float64,1}, ::OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}) at ./abstractarray.jl:798

showarg error for 0-dimensional array

julia> a = OffsetArray(reshape([1]))
0-dimensional OffsetArray(::Array{Int64,0}, Error showing value of type OffsetArray{Int64,0,Array{Int64,0}}:
ERROR: MethodError: no method matching printindices(::IOContext{Base.Terminals.TTYTerminal})
Closest candidates are:
  printindices(::IO, ::Any) at /home/fredrik/.julia/v0.7/OffsetArrays/src/OffsetArrays.jl:242
  printindices(::IO, ::Any, ::Any...) at /home/fredrik/.julia/v0.7/OffsetArrays/src/OffsetArrays.jl:240
Stacktrace:
 [1] showarg(::IOContext{Base.Terminals.TTYTerminal}, ::OffsetArray{Int64,0,Array{Int64,0}}, ::Bool) at /home/fredrik/.julia/v0.7/OffsetArrays/src/OffsetArrays.jl:236
 [2] #showarray#319(::Bool, ::Function, ::IOContext{Base.Terminals.TTYTerminal}, ::OffsetArray{Int64,0,Array{Int64,0}}, ::Bool) at ./show.jl:2063
 [3] show(::IOContext{Base.Terminals.TTYTerminal}, ::MIME{Symbol("text/plain")}, ::OffsetArray{Int64,0,Array{Int64,0}}) at ./replutil.jl:139
 [4] display(::Base.REPL.REPLDisplay{Base.REPL.LineEditREPL}, ::MIME{Symbol("text/plain")}, ::OffsetArray{Int64,0,Array{Int64,0}}) at ./repl/REPL.jl:125
 [5] display(::Base.REPL.REPLDisplay{Base.REPL.LineEditREPL}, ::OffsetArray{Int64,0,Array{Int64,0}}) at ./repl/REPL.jl:128
 [6] display(::OffsetArray{Int64,0,Array{Int64,0}}) at ./multimedia.jl:283
 [7] (::getfield(Base, Symbol("#inner#4")){NamedTuple{(),Tuple{}},typeof(display),Tuple{OffsetArray{Int64,0,Array{Int64,0}}}})() at ./essentials.jl:649
 [8] print_response(::Base.Terminals.TTYTerminal, ::Any, ::Void, ::Bool, ::Bool, ::Void) at ./repl/REPL.jl:146
 [9] print_response(::Base.REPL.LineEditREPL, ::Any, ::Void, ::Bool, ::Bool) at ./repl/REPL.jl:132
 [10] (::getfield(Base.REPL, Symbol("#do_respond#17")){Bool,getfield(Base.REPL, Symbol("##27#37")){Base.REPL.LineEditREPL,Base.REPL.REPLHistoryProvider},Base.REPL.LineEditREPL,Base.LineEdit.Prompt})(::Base.LineEdit.MIState, ::Base.GenericIOBuffer{Array{UInt8,1}}, ::Bool) at ./repl/REPL.jl:676
 [11] top-level scope

Error is here:

OffsetArrays.jl/src/OffsetArrays.jl

Line 232 in 9058120

printindices(io, indices(a)...)

because indices(a) === ().

Test failure for `@test isa(B, OffsetArray{Float32,2})`

It looks like since 0.10.0,

OffsetArrays.jl/test/runtests.jl

Line 186 in 6fa3b00

@test isa(B, OffsetArray{Float32,2})

is triggering a test failure:

$ JULIA_LOAD_PATH=src:$JULIA_LOAD_PATH julia test/runtests.jl
Skipping Base.active_repl
Skipping Base.active_repl_backend
Test Failed at /tmp/OffsetArrays.jl-0.10.0/test/runtests.jl:186
  Expression: B isa OffsetArray{Float32, 2}
   Evaluated: Float32[1.4013e-45 0.0155322; 0.0 4.56571e-41] isa OffsetArray{Float32,2,AA} where AA<:AbstractArray
ERROR: LoadError: There was an error during testing
in expression starting at /tmp/OffsetArrays.jl-0.10.0/test/runtests.jl:186

Let me know what other info you need or how I can help debug etc.

---

julia> versioninfo()
Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: AMD Ryzen 7 1700 Eight-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, znver1)

Isn't it better for OffsetArrays to have an additional method for OffsetArray argument?

Currently, when a = [1,2,3] and a_o = OffsetArray(a, -1), a_o1 = OffsetArray(a_o, 1) is of type OffsetArrays.OffsetArray{Int64,1,OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}}}. Likewise, a_o2 = OffsetArray(a_o1, -1) is of type OffsetArrays.OffsetArray{Int64,1,OffsetArrays.OffsetArray{Int64,1,OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}}}}. (Julia version: 0.6.2.2)

This hampers operability of offset arrays. Wouldn't it be better if a_o1 has the same type as a_o?

using OffsetArrays
a = [1, 2, 3]
a_o = OffsetArray(a, -1)
a_o1 = OffsetArray(a_o, 1)
a_o2 = OffsetArray(a_o1, -1)
print(typeof(a), "\n", typeof(a_o), "\n", typeof(a_o1), "\n", typeof(a_o2), "\n")

prints

Array{Int64,1}
OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}}
OffsetArrays.OffsetArray{Int64,1,OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}}}
OffsetArrays.OffsetArray{Int64,1,OffsetArrays.OffsetArray{Int64,1,OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}}}}

Is there any way to push!/pop!/delete! from an OffsetArray?

Maybe I haven't understood the usage of OffsetArray fully. I am porting some code from C++ and would like to have an 0-based indexing array on which one can still perform actions such as push! or delete! (i.e. the length is not fixed). Is this possible to achieve with this package? Or should I consider writing my own type extending AbstractArray.

The algorithm (which involves the depth of a tree) will be all over the place without 0-based indexing, which is quite annoying. Basically what I want is an array similar to the default array, with the only difference being that its index starts from 0.

Type instability in axes(arr,dim)

The method returns an IdOffsetRange or a Base.OneTo depending on whether the value of dim is less than ndims(arr). For example:

julia> a = OffsetVector(rand(3),-1:1);

julia> @code_warntype axes(a,1)
Variables
  #self#::Core.Compiler.Const(axes, false)
  A::OffsetArray{Float64,1,Array{Float64,1}}
  d::Int64

Body::Union{Base.OneTo{Int64}, OffsetArrays.IdOffsetRange{Int64,Base.OneTo{Int64}}}
1 ─       nothing
│   %2  = OffsetArrays.ndims(A)::Core.Compiler.Const(1, false)
│   %3  = (d <= %2)::Bool
└──       goto #3 if not %3
2 ─ %5  = OffsetArrays.parent(A)::Array{Float64,1}
│   %6  = OffsetArrays.axes(%5, d)::Base.OneTo{Int64}
│   %7  = Base.getproperty(A, :offsets)::Tuple{Int64}
│   %8  = Base.getindex(%7, d)::Int64
│   %9  = OffsetArrays.IdOffsetRange(%6, %8)::OffsetArrays.IdOffsetRange{Int64,Base.OneTo{Int64}}
└──       return %9
3 ─ %11 = Base.OneTo::Core.Compiler.Const(Base.OneTo, false)
│   %12 = (%11)(1)::Core.Compiler.Const(Base.OneTo(1), false)
└──       return %12

Perhaps this might return an IdOffsetRange in each case?

view != getindex, with IdentityUnitRange

I got some strange results copying a piece into a larger array, and if I understand right it looks like a problem with how view deals with IdentityUnitRange.

julia> out = reshape(1:15, 5,3) .+ 0; 

julia> item1 = 101:100:301;

julia> out[axes(item1)..., 1] .= item1; 

julia> out 
5×3 Array{Int64,2}:
 101   6  11
 201   7  12
 301   8  13
   4   9  14
   5  10  15

julia> item = OffsetArray(102:100:302, +2)
102:100:302 with indices 3:5

julia> out[axes(item)..., 1] # looks correct  
3-element OffsetArray(::Array{Int64,1}, 3:5) with eltype Int64 with indices 3:5:
 301
   4
   5

julia> view(out, axes(item)..., 1) # wrong? 
3-element view(::Array{Int64,2}, :, 1) with eltype Int64 with indices 3:5:
 5
 6
 7

julia> out[axes(item)..., 2] .= item;

julia> copyto!(view(out, axes(item)..., 2), item); # same effect

julia> out # values of item are in the wrong places!
5×3 Array{Int64,2}:
 101    6  202
 201    7  302
 301    8   13
   4    9   14
   5  102   15

(@v1.5) pkg> st OffsetArrays
Status `~/.julia/environments/v1.5/Project.toml`
  [6fe1bfb0] OffsetArrays v0.11.4

Reduction over OffsetArray error

As discussed on Discourse:
https://discourse.julialang.org/t/reducing-sum-of-offsetarray/36571

using OffsetArrays
A = rand(3,3)
B = OffsetArray(A, -1:1, 5:7)
sum(A, dims=1)   # OK
sum(B, dims=1)   # errors

Here’s the output:

Julia-1.4.0> sum(B, dims=1)
ERROR: ArgumentError: No method is implemented for reducing index range of type typeof(i).
Please implement reduced_index for this index type or report this as an issue.

Stacktrace:
 [1] reduced_index(::OffsetArrays.IdOffsetRange{Int64,Base.OneTo{Int64}}) at .\reducedim.jl:8
 [2] reduced_indices(::Tuple{OffsetArrays.IdOffsetRange{Int64,Base.OneTo{Int64}},OffsetArrays.IdOffsetRange{Int64,Base.OneTo{Int64}}}, ::Int64) at .\reducedim.jl:23
 [3] reduced_indices at .\reducedim.jl:15 [inlined]
 [4] reducedim_initarray(::OffsetArray{Float64,2,Array{Float64,2}}, ::Int64, ::Float64, ::Type{Float64}) at .\reducedim.jl:92
 [5] reducedim_initarray at .\reducedim.jl:93 [inlined]
 [6] reducedim_init at .\reducedim.jl:172 [inlined]
 [7] _mapreduce_dim at .\reducedim.jl:317 [inlined]
 [8] #mapreduce#580 at .\reducedim.jl:307 [inlined]
 [9] _sum at .\reducedim.jl:679 [inlined]
 [10] _sum at .\reducedim.jl:678 [inlined]
 [11] #sum#583 at .\reducedim.jl:652 [inlined]
 [12] top-level scope at REPL[5]:1

Suggested fix by @AndiMD

Base.reduced_index(i::OffsetArrays.IdOffsetRange) = oftype(i, first(i):first(i))

Broadcast with 1-tuple broke between 1.0.2 and 1.0.3

On Julia v1.0.2:

julia> using OffsetArrays
julia> OffsetArray([1, 2], -1) .+ (1,)
OffsetArray(::Array{Int64,1}, 0:1) with eltype Int64 with indices 0:1:
 2
 3

On Julia v1.0.3:

julia> using OffsetArrays
julia> OffsetArray([1, 2], -1) .+ (1,)
ERROR: ArgumentError: first element must be 1, got 0
Stacktrace:
 [1] (::getfield(Base, Symbol("#throwstart#30")))(::Base.Slice{UnitRange{Int64}}) at ./range.jl:291
 [2] _bcs1 at ./range.jl:293 [inlined]
 [3] _bcs at ./broadcast.jl:433 [inlined]
 [4] broadcast_shape at ./broadcast.jl:427 [inlined]
 [5] combine_axes at ./broadcast.jl:422 [inlined]
 [6] instantiate at ./broadcast.jl:266 [inlined]
 [7] materialize(::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(+),Tuple{OffsetArray{Int64,1,Array{Int64,1}},Tuple{Int64}}}) at ./broadcast.jl:756
 [8] top-level scope at none:0