Fortran-ähnliche Arrays wie FArray (Float64, -1: 1, -7: 7, -128: 512) in Julia

Im Allgemeinen ist es eine gute Entscheidung, ein 1-basiertes Array für Julia zu haben, aber manchmal ist es wünschenswert, ein Fortran-ähnliches Array mit Indizes zu haben, die einige Unterbereiche von ℤ umfassen:

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

wo es sinnvoll wäre:

in den Codes zum Buch Numerische Lösung hyperbolischer partieller Differentialgleichungen von prof. John A. Trangenstein Diese negativen Indizes werden intensiv für Geisterzellen für Randbedingungen verwendet. Gleiches gilt für Clawpack (steht für „Conservation Laws Package“) von prof. Randall J. LeVequehttp://depts.washington.edu/clawpack/und es gibt viele andere Codes, in denen solche Indizes natürlich wären. Eine solche Hilfsklasse wäre also nützlich, um solche Codes schnell zu übersetzen.

Ich habe gerade damit begonnen, einen solchen Hilfstyp zu implementieren, aber da ich für Julia noch ziemlich neu bin, wäre Ihre Hilfe sehr dankbar.

Ich habe angefangen mit:

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

Ausgabe:

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] =
...

Da ich für Julia völlig neu bin, wären alle Empfehlungen, die zu mehr Effizienz führen, sehr dankbar.

[Bearbeiten]

Das Thema wurde hier diskutiert:

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

Während der Diskussion wurden zwei Möglichkeiten ausgearbeitet, Julia-Arrays mit beliebiger Basis zu haben: SubArray-basierte Beispielnutzung mit einer Hilfsfunktion:

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

Bitte beachten Sie, dass die Grenzen nicht vollständig überprüft werden. Weitere Details finden Sie hier:https://github.com/JuliaLang/julia/issues/4044

Im Moment hat SubArray Leistungsprobleme, aber möglicherweise wird seine Leistung verbessert, siehe auch:

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

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

Ein weiterer Ansatz, der im Moment eine bessere Leistung aufweist, überprüft außerdem beide Grenzen:

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 und setindex! Methoden für FArray wurden vom Code base / array.jl inspiriert.

Anwendungsfälle:

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