AnyMath 文档

亚阵列;亚阵列

朱莉娅的 子阵列 类型是编码父级"视图"的容器 抽象阵列. 本页记录了一些设计原则和实现 子阵列s.

其中一个主要设计目标是确保两个视图的高性能 索引线[医索引]数组。 此外, 索引线 数组本身应该是 索引线 到可能的程度。

索引替换

考虑制作3d阵列的2d切片:

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

查看 删除"单例"维度(由 Int型),所以两者都 中一中二 是二维的 子阵列因此,索引这些的自然方法是 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       # for linear indexing and pointer, only valid when L==true
    stride1::Int       # used only for linear indexing
    ...
end

子阵列 有5个类型参数。 前两个是标准元素类型和维度。 接下来是父级的类型 抽象阵列. 最常用的是第四个参数,a 元组 每个维度的索引类型。 最后一个, L,只是为了方便调度而提供的;它是一个布尔值,表示索引类型是否支持快速线性索引。 稍后再谈。

如果在我们上面的例子中 A 是一个 阵列{Float64, 3},我们的 中一 上述情况将是 子阵列{Float64,2,Array{Float64,3},元组{Base.Slice{Base.OneTo{Int64}},Int64,单位范围{Int64}},假}. 请特别注意tuple参数,它存储用于创建索引的类型 中一. 同样,

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

存储这些值允许索引替换,并且将类型编码为参数允许分派到有效的算法。

索引翻译

执行索引转换需要您为不同的对象执行不同的操作 子阵列 类型。 例如,对于 中一,需要应用 我,j 索引到父数组的第一个和第三个维度,而对于 中二 人们需要将它们应用于第二和第三。 最简单的索引方法是在运行时进行类型分析:

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

不幸的是,这在性能方面将是灾难性的:每个元素访问都会分配内存,并涉及大量类型不佳的代码的运行。

更好的方法是分派到特定的方法来处理每种类型的存储索引。 就是这样 重新索引 does:它调度第一个存储索引的类型,并消耗适当数量的输入索引,然后对其余索引进行递归。 的情况下 中一,这扩展到

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

对于任何一对指数 (i,j) (除 [医CartesianIndex]s及其阵列,见下文)。

这是一个核心 子阵列 索引方法取决于 重新索引 来做这个索引翻译。 然而,有时候,我们可以避免间接,并使其更快。

线性索引

线性索引可以有效地实现时,整个数组具有分离连续元素的单个步幅,从一些偏移开始。 这意味着我们可以预先计算这些值,并将线性索引简单地表示为加法和乘法,避免了 重新索引 而且(更重要的是)笛卡尔坐标的计算速度很慢。

子阵列 类型,高效线性索引的可用性纯粹基于索引的类型,而不依赖于像父数组大小这样的值。 你可以问一组给定的索引是否支持快速线性索引与内部 基地。查看索引 功能:

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

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

这是在建造 子阵列 并存储在 L 将参数类型为布尔值,用于编码快速线性索引支持。 虽然不是绝对必要的,但这意味着我们可以直接在 子阵列{T,N,A,I,true} 没有任何中介。

由于此计算不依赖于运行时值,因此可能会错过一些步幅恰好是一致的情况:

julia> A = reshape(1:4&ast;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

构造为 视图(A,2:2:4,:) 恰好有均匀的步幅,因此线性索引确实可以有效地执行。 但是,在这种情况下,成功取决于数组的大小:如果第一个维度是奇数,

julia> A = reshape(1:5&ast;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,:] 没有统一的步幅,因此我们不能保证高效的线性索引。 因为我们必须完全基于在参数中编码的类型来做出这个决定。 子阵列, S=视图(A,2:2:4,:) 无法实现高效的线性索引。

一些细节

*请注意, 基地。重新索引 函数与输入索引的类型无关;它只是决定如何以及在哪里重新索引存储的索引。 它不仅支持整数索引,还支持非标量索引。 这意味着视图的视图不需要两个级别的间接;它们可以简单地将索引重新计算到原始父数组中! *希望现在已经相当清楚了,支持切片意味着由参数给出的维数 N,不一定等于父数组的维数或 指数 元组。 用户提供的索引也不一定与 指数 元组(例如,第二个用户提供的索引可能对应于父数组的第三维,以及 指数 元组)。 可能不太明显的是,存储的父数组的维数必须等于 指数 元组。 一些例子:

+

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.家长=AS.指数= (:,:,1:1),但支持这会使重新索引过程显着复杂化,特别是对于视图的视图。 您不仅需要对存储索引的类型进行调度,还需要检查给定索引是否是最终索引,并将任何剩余的存储索引"合并"在一起。 这不是一件容易的事,更糟糕的是:它很慢,因为它隐含地依赖于线性索引。 幸运的是,这正是计算 重新映射阵列 执行,如果可能的话,它是线性的。 因此, 查看 通过在需要时对父数组进行重新整形,确保父数组是给定索引的适当维数。 内在 子阵列 构造函数确保满足此不变量。

* [医CartesianIndex]和阵列扔一个讨厌的扳手到 重新索引 计划。 回想一下 重新索引 只需调度存储的索引的类型,以确定应该使用多少传递的索引以及它们应该去哪里。 但随着 [医]CartesianIndex,传递的参数的数量和它们索引到的维度的数量之间不再是一对一的对应关系。 如果我们回到上面的例子 基地。reindex(S1,S1。指数(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)])

+ 正确执行此操作将需要_combined_dispatch以难以处理的方式跨所有维度组合存储和传递索引。 因此, 重新索引 决不能被召唤 [医]CartesianIndex 指数。 幸运的是,标量情况很容易通过首先扁平化 [医]CartesianIndex 普通整数的参数。 数组的 [医]CartesianIndex 但是,不能如此容易地分割成正交的部分。 在尝试使用之前 重新索引, 查看 必须确保没有 [医]CartesianIndex 在参数列表中。 如果有的话,它可以简单地通过避免"平底船" 重新索引 计算完全,构造一个嵌套 子阵列 用两个层次的间接代替。