Fortran-like arrays such as FArray(Float64, -1:1,-7:7,-128:512) in Julia

Question

Generally having 1-based array for Julia is a good decision, but sometimes it is desirable to have Fortran-like array with indices that span some subranges of ℤ:

julia> x = FArray(Float64, -1:1,-7:7,-128:512)

where it would be useful:

in the codes accompanying the book Numerical Solution of Hyperbolic Partial Differential Equations by prof. John A. Trangenstein these negative indices are used intensively for ghost cells for boundary conditions. The same is true for Clawpack (stands for “Conservation Laws Package”) by prof. Randall J. LeVeque http://depts.washington.edu/clawpack/ and there are many other codes where such indices would be natural. So such auxiliary class would be useful for speedy translation of such codes.

I just started to implement such an auxiliary type but as I'm quite new to Julia your help would be greatly appreciated.

I started with:

type FArray
    ranges
    array::Array
    function FArray(T, r::Range1{Int}...)
        dims = map((x) -> length(x), r)
        array = Array(T, dims)
        new(r, array)
    end
end

Output:

julia> include ("FortranArray.jl")
julia> x = FArray(Float64, -1:1,-7:7,-128:512)
FArray((-1:1,-7:7,-128:512),3x15x641 Array{Float64,3}:
[:, :, 1] =
6.90321e-310  2.6821e-316   1.96042e-316  0.0  0.0  0.0  9.84474e-317  …  1.83233e-316  2.63285e-316  0.0           9.61618e-317  0.0        
6.90321e-310  6.32404e-322  2.63285e-316  0.0  0.0  0.0  2.63292e-316     2.67975e-316
...
[:, :, 2] =
...

As I'm completely new to Julia any recommendations especially that lead to more efficient would be greatly appreciated.

[Edit]

The topic has been discussed here:

https://groups.google.com/forum/#!topic/julia-dev/NOF6MA6tb9Y

During the discussion two ways to have Julia arrays with arbitrary base were elaborated: SubArray-based, sample usage is with a helper function:

function farray(T, r::Range1{Int64}...)
    dims = map((x) -> length(x), r)
    array = Array(T, dims)
    sub_indices = map((x) -> -minimum(x) + 2 : maximum(x), r)
    sub(array, sub_indices)
end

julia> y[-1,-7,-128] = 777
777

julia> y[-1,-7,-128] + 33
810.0

julia> y[-2,-7,-128]
ERROR: BoundsError()
 in getindex at subarray.jl:200

julia> y[2,-7,-128]
2.3977385e-316

Please note, that bounds are not checked fully more details are here: https://github.com/JuliaLang/julia/issues/4044

At the moment SubArray has performance issues but eventually its performance might be improved, see also:

https://github.com/JuliaLang/julia/issues/5117

https://github.com/JuliaLang/julia/issues/3496

Another approach that has better performance at the moment, besides checks both bounds:

type FArray{T<:Number, N, A<:AbstractArray} <: AbstractArray

    ranges
    offsets::NTuple{N,Int}
    array::A

    function FArray(r::Range1{Int}...)
        dims = map((x) -> length(x), r)
        array = Array(T, dims)
        offs = map((x) -> 1 - minimum(x), r)
        new(r, offs, array)
    end
end

FArray(T, r::Range1{Int}...) = FArray{T, length(r,), Array{T, length(r,)}}(r...)

getindex{T<:Number}(FA::FArray{T}, i1::Int) = FA.array[i1+FA.offsets[1]]
getindex{T<:Number}(FA::FArray{T}, i1::Int, i2::Int) = FA.array[i1+FA.offsets[1], i2+FA.offsets[2]]
getindex{T<:Number}(FA::FArray{T}, i1::Int, i2::Int, i3::Int) = FA.array[i1+FA.offsets[1], i2+FA.offsets[2], i3+FA.offsets[3]]

setindex!{T}(FA::FArray{T}, x, i1::Int) = arrayset(FA.array, convert(T,x), i1+FA.offsets[1])
setindex!{T}(FA::FArray{T}, x, i1::Int, i2::Int) = arrayset(FA.array, convert(T,x), i1+FA.offsets[1], i2+FA.offsets[2])
setindex!{T}(FA::FArray{T}, x, i1::Int, i2::Int, i3::Int) = arrayset(FA.array, convert(T,x), i1+FA.offsets[1], i2+FA.offsets[2], i3+FA.offsets[3])

getindex and setindex! methods for FArray were inspired by base/array.jl code.

Use cases:

julia> y = FArray(Float64, -1:1,-7:7,-128:512);

julia> y[-1,-7,-128] = 777
777

julia> y[-1,-7,-128] + 33
810.0

julia> y[-1,2,3]
0.0

julia> y[-2,-7,-128]
ERROR: BoundsError()
 in getindex at FortranArray.jl:27

julia> y[2,-7,-128]
ERROR: BoundsError()
 in getindex at FortranArray.jl:27

Answer 1

There are now two packages that provide this kind of functionality. For arrays with arbitrary start indices, see https://github.com/alsam/OffsetArrays.jl. For even more flexibility see https://github.com/mbauman/AxisArrays.jl, where indices do not have to be integers.