SubArrays

Тип SubArray в Julia является контейнером, который кодирует "представление" родительского AbstractArray. Эта страница документирует некоторые из принципов проектирования и реализации SubArray.

Одной из основных целей проектирования является обеспечение высокой производительности для представлений как IndexLinear, так и IndexCartesian массивов. Более того, представления массивов IndexLinear должны сами быть IndexLinear в той мере, в какой это возможно.

Index replacement

Рассмотрите возможность создания 2d срезов 3d массива:

julia> A = rand(2,3,4);

julia> S1 = view(A, :, 1, 2:3)
2×2 view(::Array{Float64, 3}, :, 1, 2:3) with eltype Float64:
 0.839622  0.711389
 0.967143  0.103929

julia> S2 = view(A, 1, :, 2:3)
3×2 view(::Array{Float64, 3}, 1, :, 2:3) with eltype Float64:
 0.839622  0.711389
 0.789764  0.806704
 0.566704  0.962715

view удаляет "одиночные" размеры (те, которые указаны с помощью Int), поэтому как S1, так и S2 являются двумерными SubArray. Следовательно, естественный способ индексирования этих массивов — это S1[i,j]. Чтобы извлечь значение из родительского массива A, естественный подход — заменить S1[i,j] на A[i,1,(2:3)[j]], а S2[i,j] на A[1,i,(2:3)[j]].

Ключевой особенностью дизайна SubArrays является то, что эта замена индекса может быть выполнена без каких-либо накладных расходов во время выполнения.

SubArray design

Type parameters and fields

Стратегия, принятая в первую очередь, выражается в определении типа:

struct SubArray{T,N,P,I,L} <: AbstractArray{T,N}
    parent::P
    indices::I
    offset1::Int       # for linear indexing and pointer, only valid when L==true
    stride1::Int       # used only for linear indexing
    ...
end

SubArray имеет 5 параметров типа. Первые два - это стандартный тип элемента и размерность. Следующий - это тип родительского AbstractArray. Наиболее часто используемый - это четвертый параметр, Tuple типов индексов для каждого измерения. Последний, L, предоставляется только для удобства диспетчеризации; это логическое значение, которое представляет, поддерживают ли типы индексов быструю линейную индексацию. Об этом позже.

Если в нашем примере выше A является Array{Float64, 3}, наш случай S1 выше будет SubArray{Float64,2,Array{Float64,3},Tuple{Base.Slice{Base.OneTo{Int64}},Int64,UnitRange{Int64}},false}. Обратите внимание на параметр кортежа, который хранит типы индексов, использованных для создания S1. Точно так же,

julia> S1.indices
(Base.Slice(Base.OneTo(2)), 1, 2:3)

Хранение этих значений позволяет замену индексов, а наличие типов, закодированных в качестве параметров, позволяет вызывать эффективные алгоритмы.

Index translation

Выполнение перевода индексов требует выполнения различных действий для различных конкретных типов SubArray. Например, для S1 необходимо применить индексы i,j к первому и третьему измерениям родительского массива, тогда как для S2 их нужно применить ко второму и третьему. Самый простой подход к индексации заключается в том, чтобы проводить анализ типа во время выполнения:

parentindices = Vector{Any}()
for thisindex in S.indices
    ...
    if isa(thisindex, Int)
        # Don't consume one of the input indices
        push!(parentindices, thisindex)
    elseif isa(thisindex, AbstractVector)
        # Consume an input index
        push!(parentindices, thisindex[inputindex[j]])
        j += 1
    elseif isa(thisindex, AbstractMatrix)
        # Consume two input indices
        push!(parentindices, thisindex[inputindex[j], inputindex[j+1]])
        j += 2
    elseif ...
end
S.parent[parentindices...]

К сожалению, это было бы катастрофично с точки зрения производительности: каждый доступ к элементу выделял бы память и включал бы выполнение большого количества плохо типизированного кода.

Лучший подход заключается в том, чтобы отправлять запросы к конкретным методам для обработки каждого типа хранимого индекса. Именно это и делает reindex: он отправляет запрос в зависимости от типа первого хранимого индекса и использует соответствующее количество входных индексов, а затем рекурсивно обрабатывает оставшиеся индексы. В случае S1 это расширяется до

Base.reindex(S1, S1.indices, (i, j)) == (i, S1.indices[2], S1.indices[3][j])

для любой пары индексов (i,j) (за исключением CartesianIndex и массивов thereof, см. ниже).

Это основа SubArray; методы индексации зависят от reindex, чтобы выполнить эту трансляцию индекса. Однако иногда мы можем избежать косвенной ссылки и сделать это еще быстрее.

Linear indexing

Линейная индексация может быть эффективно реализована, когда весь массив имеет одну единственную шагу, которая отделяет последовательные элементы, начиная с некоторого смещения. Это означает, что мы можем заранее вычислить эти значения и представить линейную индексацию просто как сложение и умножение, избегая косвенной индексации reindex и (что более важно) медленного вычисления декартовых координат полностью.

Для типов SubArray доступность эффективной линейной индексации основана исключительно на типах индексов и не зависит от значений, таких как размер родительского массива. Вы можете узнать, поддерживает ли данный набор индексов быструю линейную индексацию, с помощью внутренней функции Base.viewindexing:

julia> Base.viewindexing(S1.indices)
IndexCartesian()

julia> Base.viewindexing(S2.indices)
IndexLinear()

Это вычисляется во время создания SubArray и хранится в параметре типа L в виде булевого значения, которое кодирует поддержку быстрого линейного индексирования. Хотя это не строго необходимо, это означает, что мы можем определить диспетчеризацию непосредственно на SubArray{T,N,A,I,true} без каких-либо посредников.

Поскольку это вычисление не зависит от значений времени выполнения, оно может пропустить некоторые случаи, когда шаг оказывается равномерным:

julia> A = reshape(1:4*2, 4, 2)
4×2 reshape(::UnitRange{Int64}, 4, 2) with eltype Int64:
 1  5
 2  6
 3  7
 4  8

julia> diff(A[2:2:4,:][:])
3-element Vector{Int64}:
 2
 2
 2

Представление, созданное как view(A, 2:2:4, :), имеет равномерный шаг, и, следовательно, линейная индексация действительно может быть выполнена эффективно. Однако успех в этом случае зависит от размера массива: если бы первый размер был нечетным,

julia> A = reshape(1:5*2, 5, 2)
5×2 reshape(::UnitRange{Int64}, 5, 2) with eltype Int64:
 1   6
 2   7
 3   8
 4   9
 5  10

julia> diff(A[2:2:4,:][:])
3-element Vector{Int64}:
 2
 3
 2

тогда A[2:2:4,:] не имеет равномерного шага, поэтому мы не можем гарантировать эффективную линейную индексацию. Поскольку мы должны основывать это решение исключительно на типах, закодированных в параметрах SubArray, S = view(A, 2:2:4, :) не может реализовать эффективную линейную индексацию.

A few details

Обратите внимание, что функция Base.reindex не зависит от типов входных индексов; она просто определяет, как и где должны быть переиндексированы сохраненные индексы. Она поддерживает не только целочисленные индексы, но и нескалярное индексирование. Это означает, что представления представлений не нуждаются в двух уровнях косвенности; они могут просто пересчитывать индексы в исходный родительский массив!
Надеюсь, к этому моменту стало довольно ясно, что поддержка срезов означает, что размерность, заданная параметром N, не обязательно равна размерности родительского массива или длине кортежа indices. Также предоставленные пользователем индексы не обязательно совпадают с записями в кортеже indices (например, второй индекс, предоставленный пользователем, может соответствовать третьему измерению родительского массива и третьему элементу в кортеже indices).
Что может быть менее очевидным, так это то, что размерность хранимого родительского массива должна быть равна количеству эффективных индексов в кортеже indices. Вот несколько примеров:
```
A = reshape(1:35, 5, 7) # A 2d parent Array
S = view(A, 2:7)         # A 1d view created by linear indexing
S = view(A, :, :, 1:1)   # Appending extra indices is supported
```
Наивно было бы думать, что вы просто можете установить S.parent = A и S.indices = (:,:,1:1), но поддержка этого значительно усложняет процесс переиндексации, особенно для представлений представлений. Вам нужно не только отправлять по типам хранимых индексов, но и проверять, является ли данный индекс последним, и "объединять" любые оставшиеся хранимые индексы вместе. Это не простая задача, и, что еще хуже: это медленно, так как это неявно зависит от линейной индексации.
К счастью, это именно то вычисление, которое выполняет ReshapedArray, и оно делает это линейно, если это возможно. Следовательно, view гарантирует, что родительский массив имеет соответствующую размерность для заданных индексов, изменяя его форму при необходимости. Внутренний конструктор SubArray гарантирует, что это инвариант соблюдается.
CartesianIndex и массивы этого типа вносят серьезные проблемы в схему reindex. Напомню, что reindex просто распределяет по типу хранимых индексов, чтобы определить, сколько переданных индексов должно быть использовано и куда они должны быть помещены. Но с CartesianIndex больше нет однозначного соответствия между количеством переданных аргументов и количеством измерений, в которые они индексируются. Если вернуться к приведенному выше примеру Base.reindex(S1, S1.indices, (i, j)), вы можете увидеть, что расширение неверно для i, j = CartesianIndex(), CartesianIndex(2,1). Оно должно пропустить CartesianIndex() полностью и вернуть:
```
(CartesianIndex(2,1)[1], S1.indices[2], S1.indices[3][CartesianIndex(2,1)[2]])
```
Вместо этого, однако, мы получаем:
```
(CartesianIndex(), S1.indices[2], S1.indices[3][CartesianIndex(2,1)])
```
Правильное выполнение этого требует совместной диспетчеризации как по сохранённым, так и по переданным индексам во всех комбинациях размерностей неразрешимым образом. Таким образом, reindex никогда не должен вызываться с индексами CartesianIndex. К счастью, скалярный случай легко обрабатывается путём предварительного преобразования аргументов CartesianIndex в простые целые числа. Однако массивы CartesianIndex не могут быть так легко разделены на ортогональные части. Прежде чем пытаться использовать reindex, view должен убедиться, что в списке аргументов нет массивов CartesianIndex. Если такие есть, он может просто "отказаться" от вычисления reindex, создав вложенный SubArray с двумя уровнями косвенности вместо этого.