Подмассивы (SubArrays)

Тип SubArray Julia — это контейнер, кодирующий представление родительского типа AbstractArray. На этой странице описаны некоторые принципы проектирования и реализации подмассивов (SubArray).

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

Замена индекса

Рассмотрим возможность создания двумерных срезов трехмерного массива.

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

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

Проектирование подмассива

Параметры типов и поля

Принятая стратегия выражается, прежде всего, в определении типа:

struct SubArray{T,N,P,I,L} <: AbstractArray{T,N}
    parent::P
    indices::I
    offset1::Int       # для линейного индексирования и указателя, допустимо только при L==true
    stride1::Int       # используется только для линейного индексирования
    ...
end

Подмассив (SubArray) имеет пять параметров типа. Первые два — это стандартный тип и размерность элемента. Следующим является тип родительского 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)

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

Преобразование индексов

Для преобразования индексов требуется выполнять разные действия для разных конкретных типов подмассива (SubArray). Например, для S1 нужно применить индексы i,j к первому и третьему измерениям родительского массива, а для S2 их нужно применить ко второму и третьему. Проще всего выполнить индексирование можно было бы путем анализа типов во время выполнения.

parentindices = Vector{Any}()
for thisindex in S.indices
    ...
    if isa(thisindex, Int)
        # Не использовать ни один из входных индексов
        push!(parentindices, thisindex)
    elseif isa(thisindex, AbstractVector)
        # Использовать входной индекс
        push!(parentindices, thisindex[inputindex[j]])
        j += 1
    elseif isa(thisindex, AbstractMatrix)
        # Использовать два входных индекса
        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) и их массивов, см. ниже).

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

Линейное индексирование

Линейное индексирование может быть эффективно реализовано в случае, когда весь массив имеет один шаг, который разделяет последовательные элементы, начиная с некоторого смещения. Это означает, что мы можем предварительно вычислить эти значения и представить линейное индексирование просто как сложение и умножение, полностью исключив косвенное действие функции 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, :) не может реализовать эффективную линейную индексацию.

Несколько деталей

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

Менее очевидным может быть тот факт, что размерность хранимого родительского массива должна быть равна количеству эффективных индексов в кортеже indices. Некоторые примеры:
```
A = reshape(1:35, 5, 7) # Двумерный родительский массив
S = view(A, 2:7)         # Одномерное представление, созданное с помощью линейного индексирования
S = view(A, :, :, 1:1)   # Поддерживается добавление дополнительных индексов
```
Наивно полагать, что можно просто задать 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) с двумя уровнями косвенности.