# Stubborn 1-based LinearIndices for the custom type with non-zero based indexing

**URL:** https://discourse.julialang.org/t/stubborn-1-based-linearindices-for-the-custom-type-with-non-zero-based-indexing/40593
**Category:** New to Julia
**Tags:** question
**Created:** [June 2, 2020, 3:20am UTC](https://discourse.julialang.org/t/stubborn-1-based-linearindices-for-the-custom-type-with-non-zero-based-indexing/40593 "2020-06-02T03:20:49Z")
**Posts on this page:** 8
**Page:** 1

<div class="post-metadata">

### Author: ![goubochlepp](https://avatars.discourse-cdn.com/v4/letter/g/ba8739/32.png) [@goubochlepp](https://discourse.julialang.org/u/goubochlepp)
#### Post date: [June 2, 2020, 3:20am UTC](https://discourse.julialang.org/t/stubborn-1-based-linearindices-for-the-custom-type-with-non-zero-based-indexing/40593/1 "2020-06-02T03:20:49Z")

</div>

Greetings and apologies if this is very basic.

I am trying to create a custom type subtyped from `AbstractArray` to dispatch on it in place of a dense array at times. I’ve been following instructions in the [manual](https://docs.julialang.org/en/v1/devdocs/offset-arrays/) for specifying arrays with custom indices. I can make dereferencing and iteration work, but can’t get `pairs()` to function properly, it stubbornly produces unit-based indices, but unlike `enumerate()`, also pairs them with the non-existant elements.

Here is my MWE

```julia
import Base: length,size,iterate,getindex,axes, has_offset_axes
struct myrange <: AbstractVector{Int}  
    LStart::Int
    LStop::Int
end
IndexStyle(::Type{myrange}) = IndexLinear()
length(r::myrange) = r.LStop-r.LStart+1
size(r::myrange) = (length(r),)
function Base.getindex(r::myrange,i::Integer)
    @boundscheck r.LStart≤ i ≤ r.LStop || throw(BoundsError(r,i))
    i
end
firstindex(r::myrange) = r.LStart
lastindex(r::myrange) = r.LStop

```

The manual says it’s better to create a custom class for the indices range, but suggests it’s possible to use `UnitRange`. This is the approach I follow as I don’t need anything fancier

```julia
@inline Base.axes(r::myrange) = (UnitRange(r.LStart,r.LStop),)
Base.has_offset_axes(::myrange) = true

```

I also overloaded `similar()` as suggested in the manual, but this does not seem to matter.

Now, this seems to work:

```julia
julia> axes(z)[1]
4:8

```

But when attempting to use the class:

```julia
julia> LinearIndices(axes(z))
5-element LinearIndices{1,Tuple{UnitRange{Int64}}}:
1
2
3
4
5

```

and consequently

```julia
julia> pairs(IndexLinear(), z)
pairs(::myrange{Int64}) with 5 entriesError showing value of type Base.Iterators.Pairs{Int64,Integer,LinearIndices{1,Tuple{UnitRange{Int64}}},myrange{Int64}}:

ERROR: BoundsError: attempt to access 5-element myrange{Int64} with indices 4:8 at index [1]

```

The post mortem reveals that `pairs()` has indeed got the wrong indices from `LinearIndices`:

```julia
julia> properties(pairs(IndexLinear(), z))
2-element Array{Tuple{Symbol,AbstractArray{T,1} where T},1}:
(:data, Integer[#undef, #undef, #undef, 4, 5])
(:itr, [1, 2, 3, 4, 5])

```

But curiously

```julia
julia> LinearIndices(axes(z)).indices
(4:8,)

```

One might think that overloading `LinearIndices` for the custom class is required, but do observe this working perfectly:

```julia
zz=OffsetArray(4:8,4:8)
LinearIndices(axes(zz))
pairs(IndexLinear(), zz)

```

`@which LinearIndices(axes(zz))` shows that both `myrange` and `OffsetArray` get dispatched to the same (default) `LinearIndices` constructor (in Base at indices.jl:444)

So my questions are:

- how is it possible that `LinearIndices` carries wrong indices despite explicitly containing the correct ones in the eponymous field? What is the difference between `axes(::myrange)` and `axes(::OffsetArray)` that causes the divergence?

- what is the minimal way to modify my class to avoid the above behaviour?

I suspect there may be hints in the paragraph below, but I am not getting conversion errors, nor does this explain the mystery of `LinearIndices(::myrange)`

> By this definition, 1-dimensional arrays always use Cartesian indexing with the array’s native indices. To help enforce this, it’s worth noting that the index conversion functions will throw an error if shape indicates a 1-dimensional array with unconventional indexing (i.e., is a `Tuple{UnitRange}` rather than a tuple of `OneTo`). For arrays with conventional indexing, these functions continue to work the same as always.

---

<div class="post-metadata">

### Author: ![oheil](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oheil/32/220745_2.png) [@oheil](https://discourse.julialang.org/u/oheil)
#### Post date: [June 2, 2020, 9:33am UTC](https://discourse.julialang.org/t/stubborn-1-based-linearindices-for-the-custom-type-with-non-zero-based-indexing/40593/2 "2020-06-02T09:33:40Z")

</div>

Sorry for not answering your questions directly.  
Do you know about

> **[GitHub - JuliaArrays/OffsetArrays.jl: Fortran-like arrays with arbitrary,...](https://github.com/JuliaArrays/OffsetArrays.jl)**
>
> Fortran-like arrays with arbitrary, zero or negative starting indices. - GitHub - JuliaArrays/OffsetArrays.jl: Fortran-like arrays with arbitrary, zero or negative starting indices.

Maybe the work has already been done.

---

<div class="post-metadata">

### Author: ![oheil](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oheil/32/220745_2.png) [@oheil](https://discourse.julialang.org/u/oheil)
#### Post date: [June 2, 2020, 9:34am UTC](https://discourse.julialang.org/t/stubborn-1-based-linearindices-for-the-custom-type-with-non-zero-based-indexing/40593/3 "2020-06-02T09:34:14Z")

</div>

And of course: WELCOME to this fine community!

---

<div class="post-metadata">

### Author: ![goubochlepp](https://avatars.discourse-cdn.com/v4/letter/g/ba8739/32.png) [@goubochlepp](https://discourse.julialang.org/u/goubochlepp)
#### Post date: [June 2, 2020, 6:49pm UTC](https://discourse.julialang.org/t/stubborn-1-based-linearindices-for-the-custom-type-with-non-zero-based-indexing/40593/4 "2020-06-02T18:49:40Z")

</div>

Thank you Oheil,

I am familiar with `OffsetArrays`, as you can see I’ve used it as a comparator to diagnose the strange behaviour of `LinearIndices`. I cannot use the class for my purposes, because I don’t want a dense array, I need a generator (subclassed from `AbstractArray`) that iterates (this is easy) and produces correctly indexed `pairs()` (tricky). Internally, `IdOffsetRange`, the class for the `OffsetArrays`’ `axes`, is quite complicated (even more so than `CustomUnitRanges.URange`) to cover lots of eventualities, and I was wondering whether I can use a much simpler `UnitRange` implementation. The manual suggests that it should be possible, but as you can see above, there are some hidden obstacles. Could it be that there are snags in Julia base (as seems to be suggested by the extended help below)?

> ```
> Construction/coercion preserves the (shifted) values of the input range, but may modify the indexes if required by the specified types. 
> 
> !!! warning
> 
> In the future, *conversion* will preserve both the values and the indices, throwing an error when this is not achievable.
> 
> ```

I did not want to go down the `CustomUnitRanges` route, because I didn’t have the time to understand the nature of the associated complications (“_The reason is that this package’s range types should be private to the module that needs them; consequently you don’t want to define a module in the global namespace._”)

But it looks like eventually I will have to, if it’s not common knowledge among Julians.

---

<div class="post-metadata">

### Author: ![goubochlepp](https://avatars.discourse-cdn.com/v4/letter/g/ba8739/32.png) [@goubochlepp](https://discourse.julialang.org/u/goubochlepp)
#### Post date: [June 2, 2020, 7:05pm UTC](https://discourse.julialang.org/t/stubborn-1-based-linearindices-for-the-custom-type-with-non-zero-based-indexing/40593/5 "2020-06-02T19:05:34Z")

</div>

Actually, `CustomUnitRanges.URange` is not of much help, at least, out of the box. `Pairs()` don’t error out, but produce wrong indexing:

```julia
julia> Iterators.Pairs(4:8, LinearIndices(URange(4,8)))
pairs(::UnitRange{Int64}) with 5 entries:
  1 => 4
  2 => 5
  3 => 6
  4 => 7
  5 => 8

```

Mystery deepens 😕

---

<div class="post-metadata">

### Author: ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [June 2, 2020, 7:23pm UTC](https://discourse.julialang.org/t/stubborn-1-based-linearindices-for-the-custom-type-with-non-zero-based-indexing/40593/6 "2020-06-02T19:23:48Z")

</div>

> [@goubochlepp](#):
>
> Internally, `IdOffsetRange` , the class for the `OffsetArrays` ’ `axes` , is quite complicated (even more so than `CustomUnitRanges.URange` ) to cover lots of eventualities, and I was wondering whether I can use a much simpler `UnitRange` implementation.

This is precisely where you’re hitting the trouble. We’ve found that one of the core properties of `axes` is that they themselves must be similarly offset in order for things to function as you want it to!

In other words, a fundamental property of an array is that:

```julia
A[axes(A)...] == A

```

That’s simple enough when everything is 1-based, but another key property we’ve found helpful is that the axes of the result is simply the concatenated axes of the indices. In other words:

```julia
axes(A[I, J, K]) = (axes(I)..., axes(J)..., axes(K)...)

```

The interesting thing is that _many other_ algorithms — even those not based on nonscalar indexing — “just work” once you’ve satisfied this.

* * *

Now, I still think you’ll have the most success just using OffsetArrays — it’s the package that represents the accumulated learnings of how to best define offset arrays that are as functional as possible.

> [@goubochlepp](#):
>
> I cannot use the class for my purposes, because I don’t want a dense array, I need a generator (subclassed from `AbstractArray` ) that iterates (this is easy) and produces correctly indexed `pairs()` (tricky)

An `OffsetArray` isn’t necessarily dense; they wrap and offset _any_ `AbstractArray` subtype:

```julia
julia> r = 3:7
3:7

julia> ro = OffsetArray(r, -4)
3:7 with indices -3:1

julia> pairs(ro)
pairs(::OffsetArray{Int64,1,UnitRange{Int64}}) with 5 entries:
  -3 => 3
  -2 => 4
  -1 => 5
  0 => 6
  1 => 7

```

If you want to define custom behaviors, you can do so by specializing `f(::OffestArray{T,N,<:AbstractRange})` — this is “piracy” if done outside of the OffsetArrays package, but I’d bet the optimizations could be accepted as PRs there. If you want to do something more robust and less pirate-y, just define your custom array (that’s 1-based) and use OffsetArrays to offset it — then extending methods with `f(::OffsetArray{T,N,<:MyOneBasedRange})` will totally be safe.

---

<div class="post-metadata">

### Author: ![goubochlepp](https://avatars.discourse-cdn.com/v4/letter/g/ba8739/32.png) [@goubochlepp](https://discourse.julialang.org/u/goubochlepp)
#### Post date: [June 8, 2020, 12:43am UTC](https://discourse.julialang.org/t/stubborn-1-based-linearindices-for-the-custom-type-with-non-zero-based-indexing/40593/7 "2020-06-08T00:43:46Z")

</div>

@mbauman, thank you for your suggestions, your intuition is right:

> We’ve found that one of the core properties of axes is that they themselves must be similarly offset in order for things to function as you want it to!

However, I would formalize it differently. I don’t think `A[axes(A)...] == A` is sufficient. In my experiments, the criterion for the` axes()` to work was

```julia
map(first, map(first∘axes, axes(z)) ) == map(first, axes(z))
map(last, map(last∘axes, axes(z)) ) == map(last, axes(z))

```

or, to make it slightly more readable in the 1-dimensional case,

`first( Base.axes1( Base.axes1(z) ) ) == first( Base.axes1(z) )`

(when `Base.axes1()` is appropriately overloaded; in other words, `axes()` should be idempotent in the sense defined above)

Please see below the two examples: one of working `myrange` implementation, in which `A[axes(A)...] ≠ A` , so it’s not necessary (though `getindex()` could be easily overloaded to make it hold), and a broken implementation in which `A[axes(A)...] == A`.

This works:

```julia
import Base: length,size,iterate,getindex,axes, has_offset_axes
struct myrange_working <: AbstractUnitRange{Int} 
    LStart::Int
    LStop::Int
end
function Base.getindex(r::myrange_working,i::Int64)
    @boundscheck r.LStart≤ i ≤ r.LStop || throw(BoundsError(r,i))
    i
end

Base.first(r::myrange_working) = r.LStart
Base.last(r::myrange_working) = r.LStop
@inline Base.axes(r::myrange_working) = (r,)
@inline axes1(r::myrange_working) = axes(r)[1]

z = myrange_working(4,8)
z[axes(z)...] # 7:11
z[axes(z)...] == z # **false**

map(first,map(first∘axes,axes(z)))== map(first,axes(z)) # **true**
map(last,map(last∘axes,axes(z)))== map(last,axes(z)) # **true**
first(Base.axes1(Base.axes1(z))) == first(Base.axes1(z)) # **true**

julia> println(LinearIndices(axes(z)))
[4, 5, 6, 7, 8]

julia> println(pairs(IndexLinear(), z))
Base.Iterators.Pairs(4 => 4,5 => 5,6 => 6,7 => 7,8 => 8)

julia> collect(z)'
1×5 LinearAlgebra.Adjoint{Int64,Array{Int64,1}}:
 4 5 6 7 8

```

This does not:

```julia
import Base: length,size,iterate,getindex,axes, has_offset_axes
struct myrange_broken <: AbstractUnitRange{Int}  
    LStart::Int
    LStop::Int
end

function Base.getindex(r::myrange_broken,i::Int64)
    @boundscheck r.LStart≤ i ≤ r.LStop || throw(BoundsError(r,i))
    i
end
function Base.getindex(r::myrange_broken, ur::UnitRange{Int64})
    @boundscheck r.LStart≤ ur.start && ur.stop ≤ r.LStop || throw(BoundsError(r,ur))
    ur
end

Base.first(r::myrange_broken) = r.LStart
Base.last(r::myrange_broken) = r.LStop
@inline Base.axes(r::myrange_broken) = (UnitRange(r.LStart,r.LStop),)
@inline axes1(r::myrange_broken) = axes(r)[1]

z = myrange_broken(4,8)
z[axes(z)...] # 4:8
z[axes(z)...] == z # **true**

map(first,map(first∘axes,axes(z)))== map(first,axes(z)) # **false**
map(last,map(last∘axes,axes(z)))== map(last,axes(z)) # **false**
first(Base.axes1(axes1(z))) == first(Base.axes1(z)) # **false**

julia> println(LinearIndices(axes(z)))
[1, 2, 3, 4, 5]

julia> println(pairs(IndexLinear(), z))
ERROR: BoundsError: attempt to access 5-element myrange_broken with indices 4:8 at index [1]

```

For my purposes `myrange_working` is sufficient, and I would prefer this minimalistic implementation to the relative complexity of `OffsetArrays`; avoiding type piracy is a bonus.

By the way, I think the manual section on _Arrays with custom indices_ may need to be updated, as it omits the crucial bits needed for the custom `AbstractUnitRange` type to work: `Base.first()` and `Base.last()`. Is this a known issue?

---

<div class="post-metadata">

### Author: ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [June 8, 2020, 4:03pm UTC](https://discourse.julialang.org/t/stubborn-1-based-linearindices-for-the-custom-type-with-non-zero-based-indexing/40593/8 "2020-06-08T16:03:46Z")

</div>

I mean, sure, you can come up with all sorts of permutations where some things work and others don’t. I was just trying to provide some motivation for the general rule we use to unify arrays both offset and non.

A minimal but fully-functional identity-offset unit range implementation is here:

> <https://github.com/JuliaLang/julia/blob/17cff0392d4615d75aace27866875243792cd5aa/base/indices.jl#L365-L393>

Note that this really isn’t a #user:first-steps sort of thing — you’re in the deep end of the internals at this point. The #user:first-steps answer is to use OffsetArrays. 🙂
