Engee 文档

子阵列(SubArrays)

Julia’SubArray’类型是编码父类型表示的容器。 'AbstractArray'。 本页描述了设计和实现子阵列("子阵列")的一些原则。

主要设计目标之一是为类型数组的表示提供高性能 '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’的类型。 最常用的是第四个参数,它是每个维度的索引类型的元组。 最后一个’L’只是为了方便调度而提供的。 这是一个布尔值,用于确定索引类型是否支持快速线性索引。 有关此的更多细节稍后。

如果在上面的例子中`A’是一个’数组{Float64, 3}`,那么在我们的例子中’S1’将具有值'SubArray{Float64,2,Array{Float64,3},元组{Base.Slice{Base.OneTo{Int64}},Int64,单位范围{Int64}},false}'。 特别要注意tuple参数,它存储用于创建’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’函数的间接影响和(更重要的是)笛卡尔坐标的缓慢计算。

对于"子数组"类型,高效线性索引的可用性仅基于索引的类型,而不依赖于父数组大小等值。 您可以使用内部`Base’函数了解给定的一组索引是否支持快速线性索引。查看索引'。

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

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

结果在子阵列(SubArray)的构造期间计算,并作为布尔值存储在`L`类型参数中,该值编码对快速线性索引的支持。 虽然不是严格强制性的,但这意味着可以直接在"子阵列"中定义调度。{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, :)` 无法实现高效的线性索引。

一些细节

  • 注意’基。reindex’函数不依赖于输入索引的类型。 它简单地定义了存储的索引应该如何以及在哪里重新索引。 它不仅支持整数索引,还支持非标量索引。 这意味着表示的表示不需要两个级别的间接操作。 他们可以简单地将索引重新计算到原始父数组。

  • 也许已经足够清楚的是,切片支持意味着`N’参数设置的维度不一定等于父数组的维度或元组`indices’的长度。 用户指定的索引也不一定与`indexes`元组中的条目匹配(例如,第二个用户指定的索引可能对应于父数组的第三维和`indexes`元组中的第三个元素)。

    不太明显的可能是存储的父数组的维度应该等于`indexes`元组中有效索引的数量这一事实。 一些例子:

    A = reshape(1:35, 5, 7) # Двумерный родительский массив
    S = view(A, 2:7)         # Одномерное представление, созданное с помощью линейного индексирования
    S = view(A, :, :, 1:1)   # Поддерживается добавление дополнительных индексов

    相信你可以简单地设置`s.parent=A`和`s.indexes是天真的= (:,:,1:1)`, 但是支持这样的表达式使得重新索引过程更加困难,特别是对于视图表示。 您不仅需要分派存储索引的类型,还需要检查给定索引是否是最终索引并合并所有剩余的存储索引。 这不是一件容易的事。 而且,它很慢,因为它隐含地依赖于线性索引。

    幸运的是,这正是ReshapedArray数组执行的计算,如果可能的话,它是线性的。 因此,视图确保父数组具有指定索引的适当维度,并在必要时更改其形状。 'SubArray’的内部构造函数强制执行此不变量。

  • 'CartesianIndex'及其数组对`reindex`方案进行了令人不快的调整。 回想一下,reindex函数只是调度存储索引的类型,以确定应该使用多少传递的索引以及在哪里使用。 但是笛卡尔索引(CartesianIndex')不会在传递的参数数量和它们被索引的维度数量之间保持一对一的对应关系。 如果我们回到上面的例子’基地。reindex(S1,S1。indexes,(i,j)),我们会看到扩展名对于`i,j=CartesianIndex(),CartesianIndex(2,1)'是不正确的。 它应该完全_skip_'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’函数之前,视图必须确保参数列表中没有`Reindex`数组。 如果有的话,它可以简单地通过构造一个嵌套的子阵列(SubArray)来完全消除`reindex`计算,而不是两个级别的间接。