Фортраноподобные массивы, такие как FArray (Float64, -1: 1, -7: 7, -128: 512) в Julia
Как правило, наличие 1-го массива для Джулии - хорошее решение, но иногда желательно иметь Fortran-подобный массив с индексами, которые охватывают некоторые поддиапазоны ℤ:
julia> x = FArray(Float64, -1:1,-7:7,-128:512)
где это было бы полезно:
в кодах, сопровождающих книгу «Численное решение гиперболических уравнений в частных производных» проф. Джон А. Трангенштейн: эти отрицательные показатели интенсивно используются для призрачных клеток в граничных условиях. То же самое относится к Clawpack (расшифровывается как «Пакет законов о сохранении») проф. Рэндалл Дж. Левекhttp://depts.washington.edu/clawpack/и есть много других кодов, где такие индексы были бы естественными. Поэтому такой вспомогательный класс был бы полезен для быстрого перевода таких кодов.
Я только начал реализовывать такой вспомогательный тип, но, поскольку я совсем новичок в Юлии, ваша помощь будет принята с благодарностью.
Я начал с:
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
Выход:
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] =
...
Поскольку я абсолютно новичок в Юлии, любые рекомендации, особенно которые приведут к повышению эффективности, будут с благодарностью приняты.
[Редактировать]
Тема обсуждалась здесь:
https://groups.google.com/forum/#!topic/julia-dev/NOF6MA6tb9Y
Во время обсуждения были разработаны два способа иметь массивы Julia с произвольной базой: на основе SubArray, пример использования с помощью вспомогательной функции:
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
Пожалуйста, обратите внимание, что границы не проверены полностью более подробно здесь:https://github.com/JuliaLang/julia/issues/4044
На данный момент у SubArray есть проблемы с производительностью, но в конечном итоге его производительность может быть улучшена, см. Также:
https://github.com/JuliaLang/julia/issues/5117
https://github.com/JuliaLang/julia/issues/3496
Другой подход, который имеет лучшую производительность на данный момент, кроме проверки обеих границ:
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 и setindex! методы для FArray были вдохновлены кодом base / array.jl.
Сценарии использования:
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