Arrays with custom indices¶

Background¶

Because unconventional indexing breaks deeply-held assumptions throughout the Julia ecosystem, early adopters running code that has not been updated are likely to experience 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!)fori=1:length(src)@inboundsdest[i]=src[i]enddestend

This code implicitly assumes that vectors are indexed from 1. Previously that was a safe assumption, so this code was fine, but (depending on what types the user passes to this function) it may no longer be safe. If this code continued to work when passed a vector with non-1 indices, it would either produce an incorrect answer or it would segfault. (If you do get segfaults, to help locate the cause try running julia with the option --check-bounds=yes.)

To ensure that such errors are caught, in Julia 0.5 both length and sizeshould throw an error when passed an array with non-1 indexing. This is designed to force users of such arrays to check the code, and inspect it for whether it needs to be generalized.

Using indices for bounds checks and loop iteration¶

indices(A) (reminiscent of size(A)) returns a tuple of AbstractUnitRange objects, specifying the range of valid indices along each dimension of A. When A has unconventional indexing, the ranges may not start at 1. If you just want the range for a particular dimension d, there is indices(A,d).

Base implements a custom range type, OneTo, where OneTo(n) means the same thing as 1:n but in a form that guarantees (via the type system) that the lower index is 1. For any new AbstractArray type, this is the default returned by indices, and it indicates that this array type uses “conventional” 1-based indexing. Note that if you don’t want to be bothered supporting arrays with non-1 indexing, you can add the following line:

@assertall(x->isa(x,Base.OneTo),indices(A))

at the top of any function.

For bounds checking, note that there are dedicated functions checkbounds and checkindex 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 if A is multi-dimensional. In “true” linear indexing, the indices always range from 1:length(A). However, this raises an ambiguity for one-dimensional arrays (a.k.a., AbstractVector): does v[i] mean linear indexing, or Cartesian indexing with the array’s native indices?

For this reason, if you want to use linear indexing in an algorithm, your best option is to get the index range by calling linearindices(A). This will return indices(A,1) if A is an AbstractVector, and the equivalent of 1:length(A) otherwise.

In a sense, one can say that 1-dimensional arrays always use Cartesian indexing. To help enforce this, it’s worth noting that sub2ind(shape,i...) and ind2sub(shape,ind) 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.

Using indices and linearindices, here is one way you could rewrite mycopy!:

function mycopy!(dest::AbstractVector,src::AbstractVector)indices(dest)==indices(src)||throw(DimensionMismatch("vectors must match"))foriinlinearindices(src)@inboundsdest[i]=src[i]enddestend

Allocating storage using generalizations of similar¶

Storage is often allocated with Array{Int}(dims) or similar(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 use similar(storagetype,shape). storagetype indicates the kind of underlying “conventional” behavior you’d like, e.g., Array{Int} or BitArray or even dims->zeros(Float32,dims) (which would allocate an all-zeros array). shape is a tuple of Integer or AbstractUnitRange values, specifying the indices that you want the result to use.

Let’s walk through a couple of explicit examples. First, if A has conventional indices, then similar(Array{Int},indices(A)) would end up calling Array{Int}(size(A)), and thus return an array. If A is an AbstractArray type with unconventional indexing, then similar(Array{Int},indices(A)) should return something that “behaves like” an Array{Int} but with a shape (including indices) that matches A. (The most obvious implementation is to allocate an Array{Int}(size(A)) and then “wrap” it in a type that shifts the indices.)

Note also that similar(Array{Int},(indices(A,2),)) would allocate an AbstractVector{Int} (i.e., 1-dimensional array) that matches the indices of the columns of A.

Deprecations¶

In generalizing Julia’s code base, at least one deprecation was unavoidable: earlier versions of Julia defined first(::Colon)=1, meaning that the first index along a dimension indexed by : is 1. This definition can no longer be justified, so it was deprecated. There is no provided replacement, because the proper replacement depends on what you are doing and might need to know more about the array. However, it appears that many uses of first(::Colon) are really about computing an index offset; when that is the case, a candidate replacement is:

indexoffset(r::AbstractVector)=first(r)-1indexoffset(::Colon)=0

In other words, while first(:) does not itself make sense, in general you can say that the offset associated with a colon-index is zero.

Do not implement size or length¶

Perhaps the majority of pre-existing code that uses size will not work properly for arrays with non-1 indices. For that reason, it is much better to avoid implementing these methods, and use the resulting MethodError to identify code that needs to be audited and perhaps generalized.

Do not annotate bounds checks¶

Julia 0.5 includes @boundscheck to annotate code that can be removed for callers that exploit @inbounds. Initially, it seems far preferable to run with bounds checking always enabled (i.e., omit the @boundscheck annotation so the check always runs).

Custom AbstractUnitRange types¶

If you’re writing a non-1 indexed array type, you will want to specialize indices so it returns a UnitRange, or (perhaps better) a custom AbstractUnitRange. The advantage of a custom type is that it “signals” the allocation type for functions like similar. If we’re writing an array type for which indexing will start at 0, we likely want to begin by creating a new AbstractUnitRange, ZeroRange, where ZeroRange(n) is equivalent to 0:n-1.

In general, you should probably not export ZeroRange from your package: there may be other packages that implement their own ZeroRange, and having multiple distinct ZeroRange types is (perhaps counterintuitively) an advantage: ModuleA.ZeroRange indicates that similar should create a ModuleA.ZeroArray, whereas ModuleB.ZeroRange indicates a ModuleB.ZeroArray type. This design allows peaceful coexistence among many different custom array types.

Note that the Julia package CustomUnitRanges.jl can sometimes be used to avoid the need to write your own ZeroRange type.

Specializing indices¶

Once you have your AbstractUnitRange type, then specialize indices:

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

where here we imagine that ZeroArray has a field called size (there would be other ways to implement this).

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

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

may not be what you want: you may need to specialize it to return something other than OneTo(1) when d>ndims(A). Likewise, in Base there is a dedicated function indices1 which is equivalent to indices(A,1) but which avoids checking (at runtime) whether ndims(A)>0. (This is purely a performance optimization.) It is defined as:

indices1{T}(A::AbstractArray{T,0})=OneTo(1)indices1{T}(A::AbstractArray{T})=indices(A)[1]

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

Specializing similar¶

Given your custom ZeroRange type, then you should also add the following two specializations for similar:

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.

Specializing reshape¶

Optionally, define a method

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

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