# Multidimensional arrays with 0-based indices

**URL:** <https://discourse.julialang.org/t/multidimensional-arrays-with-0-based-indices/3325>\
**Category:** General Usage\
**Tags:** question\
**Created:** [April 21, 2017, 12:19pm UTC](https://discourse.julialang.org/t/multidimensional-arrays-with-0-based-indices/3325 "2017-04-21T12:19:23Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![barche](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barche/32/79_2.png) [@barche](https://discourse.julialang.org/u/barche)\
**Post date:** [April 21, 2017, 12:19pm UTC](https://discourse.julialang.org/t/multidimensional-arrays-with-0-based-indices/3325/1 "2017-04-21T12:19:23Z")

</div>

Hi,

I’m trying to create an array view (for data from C++) that uses 0-based indices. I got most of it working, except for the use of `end` as an index into a multidimensional error. The following code mimics what I want to do in a standalone example, the final test fails:

```julia
using Compat
using Base.Test

using CustomUnitRanges: filename_for_zerorange
include(filename_for_zerorange)

immutable PtrView{ScalarT,N,LayoutT} <: AbstractArray{ScalarT,N}
  ptr::Ptr{ScalarT}
  size::NTuple{N,Int}
end

immutable LayoutLeft
end

PtrView{ScalarT,N,LayoutT}(A::Array{ScalarT,N},::Type{LayoutT}) = PtrView{ScalarT,N,LayoutT}(pointer(A), size(A))

@compat Base.IndexStyle{ScalarT,N,LayoutT}(::Type{PtrView{ScalarT,N,LayoutT}}) = IndexLinear()
Base.indices{ScalarT,N}(v::PtrView{ScalarT,N,LayoutLeft}) = map(ZeroRange, v.size)
Base.endof{ScalarT,N}(v::PtrView{ScalarT,N,LayoutLeft}) = prod(v.size)-1
Base.getindex{ScalarT,N}(v::PtrView{ScalarT,N,LayoutLeft}, i::Integer) = unsafe_load(v.ptr,i+1)
Base.setindex!{ScalarT,N}(v::PtrView{ScalarT,N,LayoutLeft}, value, i::Integer) = unsafe_store!(v.ptr,value,i+1)

v_arr = [1,2,3]
A_arr = [1 2 3;4 5 6]

v = PtrView(v_arr, LayoutLeft)
A = PtrView(A_arr, LayoutLeft)

@show @test v[0] == 1
@show @test v[end] == 3
@show @test A[0] == 1
@show @test A[end] == 6
@show @test A[0,end] == 3

```

The error:

```nohighlight
Error During Test
  Test threw an exception of type MethodError
  Expression: A[0,end] == 3
  MethodError: no method matching size(::PtrView{Int64,2,LayoutLeft})
  Closest candidates are:
    size{T,N}(::AbstractArray{T,N}, !Matched::Any) at abstractarray.jl:47
    size{N}(::Any, !Matched::Integer, !Matched::Integer, !Matched::Integer...) at abstractarray.jl:48
    size(!Matched::BitArray{1}) at bitarray.jl:39
    ...
   in size at ./abstractarray.jl:47 [inlined]
   in trailingsize(::PtrView{Int64,2,LayoutLeft}, ::Int64) at ./abstractarray.jl:201
ERROR: LoadError: There was an error during testing
 in record(::Base.Test.FallbackTestSet, ::Base.Test.Error) at ./test.jl:397
 in do_test(::Base.Test.Threw, ::Expr) at ./test.jl:281

```

Do I get around this by overriding `Base.trailingsize`? This seems related but unconcluded:  
[https://github.com/JuliaLang/julia/issues/16104](https://github.com/JuliaLang/julia/issues/16104)

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [April 21, 2017, 1:18pm UTC](https://discourse.julialang.org/t/multidimensional-arrays-with-0-based-indices/3325/2 "2017-04-21T13:18:50Z")

</div>

With unconventional indices, you’re much better off using `ibegin` and `iend`, see [GitHub - JuliaArrays/EndpointRanges.jl: Julia package for doing arithmetic on endpoints in array indexing](https://github.com/JuliaArrays/EndpointRanges.jl).

Presumably you’re aware that this is a “dangerous” view in that you need to be sure that `ptr` doesn’t go away. If it’s owned by C code and won’t be GCed, you should be safe, but some care is required.

---

<div class="post-metadata">

**Author:** ![barche](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barche/32/79_2.png) [@barche](https://discourse.julialang.org/u/barche)\
**Post date:** [April 21, 2017, 2:37pm UTC](https://discourse.julialang.org/t/multidimensional-arrays-with-0-based-indices/3325/3 "2017-04-21T14:37:26Z")

</div>

Thanks, EndpointRanges looks very useful. I looked a bit at how `end` is implemented, and at a first glance adding a function like this:

```julia
Base.endof(A::AbstractArray, d) = last(indices(A,d))

```

and then making end call that instead of `size` or `trailingsize` would fix this. Probably it’s not that simple?

Also, I noticed linearindices is defined as:

```julia
linearindices(A) = OneTo(_length(A))

```

This should also be overridden for multidimensional arrays with unconventional indices, right?

The memory is indeed managed by C++ in the real use case.

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [April 21, 2017, 2:41pm UTC](https://discourse.julialang.org/t/multidimensional-arrays-with-0-based-indices/3325/4 "2017-04-21T14:41:45Z")

</div>

I’d forgotten we had the `d` version of `endof`, so that seems good. OTOH, any code that does `A[..., 1:end, ...]` is probably not what you want (it leaves out the first element), so there’s a sense in which _not_ defining that `endof` method might be good defensive practice, at least in the early days.

---

<div class="post-metadata">

**Author:** ![barche](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barche/32/79_2.png) [@barche](https://discourse.julialang.org/u/barche)\
**Post date:** [April 21, 2017, 3:29pm UTC](https://discourse.julialang.org/t/multidimensional-arrays-with-0-based-indices/3325/5 "2017-04-21T15:29:25Z")

</div>

I don’t think here is a `d` version of `endof`, unless I missed it. I was contemplating a PR to add it, but I guess ranges complicate things. Is it worth attempting this anyway, or is something like EndpointRanges.jl going to replace this functionality in the “end”? Also, is there a test suite that checks all the range corner cases even with custom indices?

Regarding `linearindices`: should it be added to the docs that this needs to be overridden for multidimensional arrays using custom indices?

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [April 21, 2017, 6:17pm UTC](https://discourse.julialang.org/t/multidimensional-arrays-with-0-based-indices/3325/6 "2017-04-21T18:17:01Z")

</div>

Sorry, I wasn’t reading carefully enough. I guess the reason I forgot about the `d` version is we don’t have it (yet) 🙂.

> Is it worth attempting this anyway, or is something like EndpointRanges.jl going to replace this functionality in the “end”?

Perhaps the best approach would be to put it in the parser. That’s the main place where `end` is implemented now. We’d want to include `begin` too. My Scheme is very weak so I haven’t tried to tackle this myself yet.

> Also, is there a test suite that checks all the range corner cases even with custom indices?

Perhaps the best would be to copy the OffsetArrays tests?

> Regarding linearindices: should it be added to the docs that this needs to be overridden for multidimensional arrays using custom indices?

I don’t think it needs specialzation. The reason is that there are two methods:

```julia
linearindices(A) = (@_inline_meta; OneTo(_length(A)))
linearindices(A::AbstractVector) = (@_inline_meta; indices1(A))

```

The second one “does the right thing” for vectors. (`indices1` is a short-circuiting equivalent of `indices(A, 1)`.) For multidimensional arrays, linear indexing always starts with 1 no matter what indices the array has, so that first definition seems to cover the bases.

That said, if you discover reason to need to specialize it, I’d be interested in hearing why.

---

<div class="post-metadata">

**Author:** ![barche](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barche/32/79_2.png) [@barche](https://discourse.julialang.org/u/barche)\
**Post date:** [April 22, 2017, 8:09am UTC](https://discourse.julialang.org/t/multidimensional-arrays-with-0-based-indices/3325/7 "2017-04-22T08:09:37Z")

</div>

> [@tim.holy](#):
>
> For multidimensional arrays, linear indexing always starts with 1 no matter what indices the array has, so that first definition seems to cover the bases.

OK, I didn’t realise linear indices were always 1-based for ND-arrays. What is actually needed is to define get/setindex differently for the 1D and ND case, then:

```julia
Base.getindex{ScalarT}(v::PtrView{ScalarT,1,LayoutLeft}, i::Integer) = unsafe_load(v.ptr,i+1)
Base.setindex!{ScalarT}(v::PtrView{ScalarT,1,LayoutLeft}, value, i::Integer) = unsafe_store!(v.ptr,value,i+1)
Base.getindex{ScalarT,N}(v::PtrView{ScalarT,N,LayoutLeft}, i::Integer) = unsafe_load(v.ptr,i)
Base.setindex!{ScalarT,N}(v::PtrView{ScalarT,N,LayoutLeft}, value, i::Integer) = unsafe_store!(v.ptr,value,i)

```

This is necessary since e.g. `A[0,0]` is converted to `A[1]`. I wanted to specialize `linearindices` because I thought `A[0]` refers to the first matrix element, so indeed there is no need to specialize it.

---

<div class="post-metadata">

**Author:** ![barche](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barche/32/79_2.png) [@barche](https://discourse.julialang.org/u/barche)\
**Post date:** [April 24, 2017, 10:14am UTC](https://discourse.julialang.org/t/multidimensional-arrays-with-0-based-indices/3325/8 "2017-04-24T10:14:45Z")

</div>

I had some more fun with this, new code for a complete example with benchmarks and omitting the pointers for clarity:

```julia
using Compat
using Base.Test

using CustomUnitRanges: filename_for_zerorange
include(filename_for_zerorange)

immutable ZeroView{ScalarT,N} <: AbstractArray{ScalarT,N}
  a::Array{ScalarT,N}
end

@compat Base.IndexStyle{ScalarT,N}(::Type{ZeroView{ScalarT,N}}) = IndexLinear()
Base.indices{ScalarT,N}(v::ZeroView{ScalarT,N}) = map(ZeroRange, size(v.a))
# needed to solve type instability:
Base.indices{ScalarT,N}(v::ZeroView{ScalarT,N}, d) = ZeroRange(size(v.a,d))
Base.getindex{ScalarT}(v::ZeroView{ScalarT,1}, i::Integer) = v.a[i+1]
Base.setindex!{ScalarT}(v::ZeroView{ScalarT,1}, value, i::Integer) = (v.a[i+1] = value)
Base.getindex{ScalarT,N}(v::ZeroView{ScalarT,N}, i::Integer) = v.a[i] # linearindices are always 1-based for ND arrays
Base.setindex!{ScalarT,N}(v::ZeroView{ScalarT,N}, value, i::Integer) = (v.a[i] = value)

v = ZeroView([1,2,3])
A = ZeroView([1 2 3;4 5 6])

println("v:\n", v)
println("A:\n", A)

@test v[0] == 1
@test v[2] == 3
@test A[1] == 1
@test A[6] == 6
@test A[0,0] == 1
@test A[0,2] == 3
@test A[1,0] == 4
@test A[1,2] == 6

function benchmark_fill_vec(A)
  for i in linearindices(A)
    A[i] = i
  end
end

function benchmark_fill_mat(A)
  for i in indices(A,1)
    for j in indices(A,2)
      A[i,j] = i*(j+1)
    end
  end
end

const bench_size = 1000000
const ncols = 3
v_bench = ZeroView(zeros(bench_size))
A_bench = ZeroView(zeros(bench_size,ncols))

println("vector timings:")
@time benchmark_fill_vec(v_bench)
@time benchmark_fill_vec(v_bench)
@time benchmark_fill_vec(v_bench)

@test v_bench[0] == 0
@test v_bench[bench_size-1] == bench_size-1

println("matrix timings, linear indexing:")
@time benchmark_fill_vec(A_bench)
@time benchmark_fill_vec(A_bench)
@time benchmark_fill_vec(A_bench)

@test A_bench[0,0] == 1
@test A_bench[bench_size-1,2] == bench_size*ncols

println("matrix timings, [i,j] indexing:")
@time benchmark_fill_mat(A_bench)
@time benchmark_fill_mat(A_bench)
@time benchmark_fill_mat(A_bench)

@test A_bench[0,0] == 0
@test A_bench[bench_size-1,2] == (bench_size-1)*ncols

```

There is a problem with the `indices(A,d)` variant, since this introduces a type instability. Adding the specialization:

```julia
Base.indices{ScalarT,N}(v::ZeroView{ScalarT,N}, d) = ZeroRange(size(v.a,d))

```

fixes that, but may break other functionality such as ranges, since the definition in Base is:

```julia
function indices{T,N}(A::AbstractArray{T,N}, d)
    @_inline_meta
    d <= N ? indices(A)[d] : OneTo(1)
end

```

I believe the branch on N is the cause of the type instability.

Running with bounds checking off and the specialized `indices` d version I get these timings (on 0.6):

```nohighlight
v:
[1,2,3]
A:
[1 2 3; 4 5 6]
vector timings:
  0.008465 seconds (6.13 k allocations: 293.714 KB)
  0.000626 seconds (4 allocations: 160 bytes)
  0.000668 seconds (4 allocations: 160 bytes)
matrix timings, linear indexing:
  0.010662 seconds (4.06 k allocations: 192.623 KB)
  0.002969 seconds (4 allocations: 160 bytes)
  0.003170 seconds (4 allocations: 160 bytes)
matrix timings, [i,j] indexing:
  0.024965 seconds (10.71 k allocations: 424.688 KB)
  0.012448 seconds (4 allocations: 160 bytes)
  0.014230 seconds (4 allocations: 160 bytes)

```

The difference in timing between the [i,j] and linear indexing surprises me.

BTW, would it be useful to add this example to the docs (once everything is correct, of course)?