Arrays with custom indices

Conventionally, Julia's arrays are indexed starting at 1, whereas some other languages start numbering at 0, and yet others (e.g., Fortran) allow you to specify arbitrary starting indices. While there is much merit in picking a standard (i.e., 1 for Julia), there are some algorithms which simplify considerably if you can index outside the range1:size(A,d) (and not just0:size(A,d)-1, either). To facilitate such computations, Julia supports arrays with arbitrary indices.

The purpose of this page is to address the question, "what do I have to do to support such arrays in my own code?" First, let's address the simplest case: if you know that your code will never need to handle arrays with unconventional indexing, hopefully the answer is "nothing." Old code, on conventional arrays, should function essentially without alteration as long as it was using the exported interfaces of Julia. If you find it more convenient to just force your users to supply traditional arrays where indexing starts at one, you can add

Base.require_one_based_indexing(arrays...)

wherearrays... is a list of the array objects that you wish to check for anything that violates 1-based indexing.

Generalizing existing code

As an overview, the steps are:

• replace many uses ofsize withaxes
• replace1:length(A) witheachindex(A), or in some casesLinearIndices(A)
• replace explicit allocations likeArray{Int}(undef, size(B)) withsimilar(Array{Int}, axes(B))

These are described in more detail below.

Things to watch out for

Because unconventional indexing breaks many people's assumptions that all arrays start indexing with 1, there is always the chance that using such arrays will trigger errors. The most frustrating bugs would be incorrect results or segfaults (total crashes of Julia). For example, consider the following function:

function mycopy!(dest::AbstractVector, src::AbstractVector)    length(dest) == length(src) || throw(DimensionMismatch("vectors must match"))    # OK, now we're safe to use @inbounds, right? (not anymore!)    for i = 1:length(src)        @inbounds dest[i] = src[i]    end    destend

This code implicitly assumes that vectors are indexed from 1; ifdest starts at a different index thansrc, there is a chance that this code would trigger a segfault. (If you do get segfaults, to help locate the cause try running julia with the option--check-bounds=yes.)

Usingaxes for bounds checks and loop iteration

axes(A) (reminiscent ofsize(A)) returns a tuple ofAbstractUnitRange{<:Integer} objects, specifying the range of valid indices along each dimension ofA. WhenA has unconventional indexing, the ranges may not start at 1. If you just want the range for a particular dimensiond, there isaxes(A, d).

Base implements a custom range type,OneTo, whereOneTo(n) means the same thing as1:n but in a form that guarantees (via the type system) that the lower index is 1. For any newAbstractArray type, this is the default returned byaxes, and it indicates that this array type uses "conventional" 1-based indexing.

For bounds checking, note that there are dedicated functionscheckbounds andcheckindex which can sometimes simplify such tests.

Linear indexing (LinearIndices)

Some algorithms are most conveniently (or efficiently) written in terms of a single linear index,A[i] even ifA is multi-dimensional. Regardless of the array's native indices, linear indices always range from1:length(A). However, this raises an ambiguity for one-dimensional arrays (a.k.a.,AbstractVector): doesv[i] mean linear indexing , or Cartesian indexing with the array's native indices?

For this reason, your best option may be to iterate over the array witheachindex(A), or, if you require the indices to be sequential integers, to get the index range by callingLinearIndices(A). This will returnaxes(A, 1) if A is an AbstractVector, and the equivalent of1:length(A) otherwise.

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 aTuple{UnitRange} rather than a tuple ofOneTo). For arrays with conventional indexing, these functions continue to work the same as always.

Usingaxes andLinearIndices, here is one way you could rewritemycopy!:

function mycopy!(dest::AbstractVector, src::AbstractVector)    axes(dest) == axes(src) || throw(DimensionMismatch("vectors must match"))    for i in LinearIndices(src)        @inbounds dest[i] = src[i]    end    destend

Allocating storage using generalizations ofsimilar

Storage is often allocated withArray{Int}(undef, dims) orsimilar(A, args...). When the result needs to match the indices of some other array, this may not always suffice. The generic replacement for such patterns is to usesimilar(storagetype, shape).storagetype indicates the kind of underlying "conventional" behavior you'd like, e.g.,Array{Int} orBitArray or evendims->zeros(Float32, dims) (which would allocate an all-zeros array).shape is a tuple ofInteger orAbstractUnitRange values, specifying the indices that you want the result to use. Note that a convenient way of producing an all-zeros array that matches the indices of A is simplyzeros(A).

Let's walk through a couple of explicit examples. First, ifA has conventional indices, thensimilar(Array{Int}, axes(A)) would end up callingArray{Int}(undef, size(A)), and thus return an array. IfA is anAbstractArray type with unconventional indexing, thensimilar(Array{Int}, axes(A)) should return something that "behaves like" anArray{Int} but with a shape (including indices) that matchesA. (The most obvious implementation is to allocate anArray{Int}(undef, size(A)) and then "wrap" it in a type that shifts the indices.)

Note also thatsimilar(Array{Int}, (axes(A, 2),)) would allocate anAbstractVector{Int} (i.e., 1-dimensional array) that matches the indices of the columns ofA.

Writing custom array types with non-1 indexing

Most of the methods you'll need to define are standard for anyAbstractArray type, seeAbstract Arrays. This page focuses on the steps needed to define unconventional indexing.

CustomAbstractUnitRange types

If you're writing a non-1 indexed array type, you will want to specializeaxes so it returns aUnitRange, or (perhaps better) a customAbstractUnitRange. The advantage of a custom type is that it "signals" the allocation type for functions likesimilar. If we're writing an array type for which indexing will start at 0, we likely want to begin by creating a newAbstractUnitRange,ZeroRange, whereZeroRange(n) is equivalent to0:n-1.

In general, you should probablynot exportZeroRange from your package: there may be other packages that implement their ownZeroRange, and having multiple distinctZeroRange types is (perhaps counterintuitively) an advantage:ModuleA.ZeroRange indicates thatsimilar should create aModuleA.ZeroArray, whereasModuleB.ZeroRange indicates aModuleB.ZeroArray type. This design allows peaceful coexistence among many different custom array types.

Note that the Julia packageCustomUnitRanges.jl can sometimes be used to avoid the need to write your ownZeroRange type.

Specializingaxes

Once you have yourAbstractUnitRange type, then specializeaxes:

Base.axes(A::ZeroArray) = map(n->ZeroRange(n), A.size)

where here we imagine thatZeroArray has a field calledsize (there would be other ways to implement this).

In some cases, the fallback definition foraxes(A, d):

axes(A::AbstractArray{T,N}, d) where {T,N} = d <= N ? axes(A)[d] : OneTo(1)

may not be what you want: you may need to specialize it to return something other thanOneTo(1) whend > ndims(A). Likewise, inBase there is a dedicated functionaxes1 which is equivalent toaxes(A, 1) but which avoids checking (at runtime) whetherndims(A) > 0. (This is purely a performance optimization.) It is defined as:

axes1(A::AbstractArray{T,0}) where {T} = OneTo(1)axes1(A::AbstractArray) = axes(A)[1]

If the first of these (the zero-dimensional case) is problematic for your custom array type, be sure to specialize it appropriately.

Specializingsimilar

Given your customZeroRange type, then you should also add the following two specializations forsimilar:

function Base.similar(A::AbstractArray, T::Type, shape::Tuple{ZeroRange,Vararg{ZeroRange}})    # bodyendfunction Base.similar(f::Union{Function,DataType}, shape::Tuple{ZeroRange,Vararg{ZeroRange}})    # bodyend

Both of these should allocate your custom array type.

Specializingreshape

Optionally, define a method

Base.reshape(A::AbstractArray, shape::Tuple{ZeroRange,Vararg{ZeroRange}}) = ...

and you canreshape an array so that the result has custom indices.

For objects that mimic AbstractArray but are not subtypes

has_offset_axes depends on havingaxes defined for the objects you call it on. If there is some reason you don't have anaxes method defined for your object, consider defining a method

Base.has_offset_axes(obj::MyNon1IndexedArraylikeObject) = true

This will allow code that assumes 1-based indexing to detect a problem and throw a helpful error, rather than returning incorrect results or segfaulting julia.

Catching errors

If your new array type triggers errors in other code, one helpful debugging step can be to comment out@boundscheck in yourgetindex andsetindex! implementation. This will ensure that every element access checks bounds. Or, restart julia with--check-bounds=yes.

In some cases it may also be helpful to temporarily disablesize andlength for your new array type, since code that makes incorrect assumptions frequently uses these functions.