返回介绍

多维数组

发布于 2019-07-03 15:53:51 字数 14117 浏览 1151 评论 0 收藏 0

与大多数技术计算语言一样,Julia 提供原生的数组实现。 大多数技术计算语言非常重视其数组实现,但需要付出使用其它容器的代价。Julia 用同样的方式来处理数组。就像和其它用 Julia 写的代码一样,Julia 的数组库几乎完全是用 Julia 自身实现的,它的性能源自编译器。这样一来,用户就可以通过继承 AbstractArray 的方式来创建自定义数组类型。 实现自定义数组类型的更多详细信息,请参阅[manual section on the AbstractArray interface]

索引 n 维数组 A 的一般语法是:

X = A[I_1, I_2, ..., I_n]

其中每个 I_k 可以是标量整数,整数数组或任何其他[支持的索引类型]

在表达式 A[I_1, I_2, ..., I_n] 中,每个 I_k 可以是标量索引,标量索引数组,或者用 to_indices 转换成的表示标量索引数组的对象:

  1. 标量索引。默认情况下,这包括:
    • 非布尔的整数
    • CartesianIndex{N}s,其行为类似于跨越多个维度的 N 维整数元组(详见下文)
  2. 标量索引数组。这包括:
    • 整数向量和多维整数数组
    • [] 这样的空数组,它不选择任何元素
    • a:ca:b:c 的范围,从 ac(包括)选择连续或间隔的部分元素
    • 任何自定义标量索引数组,它是 AbstractArray 的子类型
    • CartesianIndex{N} 数组(详见下文)
  3. 一个表示标量索引数组的对象,可以通过to_indices转换为这样的对象。 默认情况下,这包括:
    • Colon() (:),表示整个维度内或整个数组中的所有索引
    • 布尔数组,选择其中值为 true 的索引对应的元素(更多细节见下文)

一些例子:

julia> A = reshape(collect(1:2:18), (3, 3))
3×3 Array{Int64,2}:
 1   7  13
 3   9  15
 5  11  17

julia> A[4]
7

julia> A[[2, 5, 8]]
3-element Array{Int64,1}:
  3
  9
 15

julia> A[[1 4; 3 8]]
2×2 Array{Int64,2}:
 1   7
 5  15

julia> A[[]]
0-element Array{Int64,1}

julia> A[1:2:5]
3-element Array{Int64,1}:
 1
 5
 9

julia> A[2, :]
3-element Array{Int64,1}:
  3
  9
 15

julia> A[:, 3]
3-element Array{Int64,1}:
 13
 15
 17

笛卡尔索引

特殊的 CartesianIndex{N} 对象表示一个标量索引,其行为类似于张成多个维度的 N 维整数元组。例如:

julia> A = reshape(1:32, 4, 4, 2);

julia> A[3, 2, 1]
7

julia> A[CartesianIndex(3, 2, 1)] == A[3, 2, 1] == 7
true

如果单独考虑,这可能看起来相对微不足道;CartesianIndex 只是将多个整数聚合成一个表示单个多维索引的对象。 但是,当与其他索引形式和迭代器组合产生多个 CartesianIndex 时,这可以生成非常优雅和高效的代码。请参阅下面的迭代,有关更高级的示例,请参阅关于多维算法和迭代博客文章

也支持 CartesianIndex {N} 的数组。它们代表一组标量索引,每个索引都跨越 N 个维度,从而实现一种有时也称为逐点索引的索引形式。例如,它可以从上面的 A 的第一「页」访问对角元素:

julia> page = A[:,:,1]
4×4 Array{Int64,2}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> page[[CartesianIndex(1,1),
             CartesianIndex(2,2),
             CartesianIndex(3,3),
             CartesianIndex(4,4)]]
4-element Array{Int64,1}:
  1
  6
 11
 16

这可以通过 [dot broadcasting](http://127.0.0.5/@ref man-vectorized) 以及普通整数索引(而不是把从 A 中提取第一“页”作为单独的步骤)更加简单地表达。它甚至可以与 : 结合使用,同时从两个页面中提取两个对角线:

julia> A[CartesianIndex.(axes(A, 1), axes(A, 2)), 1]
4-element Array{Int64,1}:
  1
  6
 11
 16

julia> A[CartesianIndex.(axes(A, 1), axes(A, 2)), :]
4×2 Array{Int64,2}:
  1  17
  6  22
 11  27
 16  32

!!! warning

`CartesianIndex` 和 `CartesianIndex` 数组与表示维度的最后一个索引的 `end` 关键字不兼容。 不要在可能包含 `CartesianIndex` 或其数组的索引表达式中使用 `end`。

逻辑索引

通常称为逻辑索引或使用逻辑掩码索引,通过布尔数组进行索引选择索引处其值为 true 的元素。 通过布尔向量 B 进行索引实际上与通过findall(B)返回的整数向量进行索引相同。 类似地,通过 N 维布尔数组进行索引实际上与通过 CartesianIndex{N} 向量在值为 true 处进行索引相同。 逻辑索引必须是与其索引的维度长度相同的向量,或者它必须是提供的唯一索引,并且与索引的数组的大小和维度相匹配。 将布尔数组直接用作索引,通常比首先调用findall更有效。

julia> x = reshape(1:16, 4, 4)
4×4 reshape(::UnitRange{Int64}, 4, 4) with eltype Int64:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> x[[false, true, true, false], :]
2×4 Array{Int64,2}:
 2  6  10  14
 3  7  11  15

julia> mask = map(ispow2, x)
4×4 Array{Bool,2}:
  true  false  false  false
  true  false  false  false
 false  false  false  false
  true   true  false   true

julia> x[mask]
5-element Array{Int64,1}:
  1
  2
  4
  8
 16

迭代

迭代整个数组的推荐方法是

for a in A
    # Do something with the element a
end

for i in eachindex(A)
    # Do something with i and/or A[i]
end

当你需要每个元素的值而不是索引时,使用第一个构造。 在第二个构造中,如果 A 是具有快速线性索引的数组类型,i 将是 Int; 否则,它将是一个 CartesianIndex

julia> A = rand(4,3);

julia> B = view(A, 1:3, 2:3);

julia> for i in eachindex(B)
           @show i
       end
i = CartesianIndex(1, 1)
i = CartesianIndex(2, 1)
i = CartesianIndex(3, 1)
i = CartesianIndex(1, 2)
i = CartesianIndex(2, 2)
i = CartesianIndex(3, 2)

for i = 1:length(A) 相比,eachindex 提供了一种迭代任何数组类型的有效方法。

Array traits

如果你编写一个自定义的 AbstractArray 类型,你可以用以下代码指定它使用快速线性索引

Base.IndexStyle(::Type{<:MyArray}) = IndexLinear()

此设置将导致 myArray 上的 eachindex 迭代使用整数。如果未指定此特征,则使用默认值 IndexCartesian()

数组、向量化操作符及函数

以下运算符支持对数组操作

  1. 一元运算符 -- -, +
  2. 二元运算符 -- -, +, *, /, \, ^
  3. 比较操作符 -- ==, !=, (isapprox),

另外,为了便于数学上和其他运算的向量化,Julia [提供了点语法(dot syntax)](http://127.0.0.5/@ref man-vectorized) f.(args...),例如,sin.(x)min.(x,y),用于数组或数组和标量的混合上的按元素运算(广播运算);当与其他点调用(dot call)结合使用时,它们的额外优点是能「融合」到单个循环中,例如,sin.(cos.(x))

此外,每个二元运算符支持相应的[点操作版本](http://127.0.0.5/@ref man-dot-operators),可以应用于此类[融合 broadcasting 操作](http://127.0.0.5/@ref man-vectorized)的数组(以及数组和标量的组合),例如 z .== sin.(x .* y)

请注意,类似 == 的比较运算在作用于整个数组时,得到一个布尔结果。使用像 .== 这样的点运算符进行按元素的比较。(对于像 < 这样的比较操作,只有按元素运算的版本 .< 适用于数组。)

还要注意 max.(a,b)maximum(a) 之间的区别,max.(a,b)ab 的每个元素 broadcasts maxmaximum(a) 寻找在 a 中的最大值。min.(a,b)minimum(a) 也有同样的关系。

广播

有时需要在不同尺寸的数组上执行元素对元素的操作,例如将矩阵的每一列加一个向量。一种低效的方法是将向量复制成矩阵的大小:

julia> a = rand(2,1); A = rand(2,3);

julia> repeat(a,1,3)+A
2×3 Array{Float64,2}:
 1.20813  1.82068  1.25387
 1.56851  1.86401  1.67846

当维度较大的时候,这种方法将会十分浪费,所以 Julia 提供了广播 broadcast,它将会将参数中低维度的参数扩展,使得其与其他维度匹配,且不会使用额外的内存,并将所给的函数逐元素地应用。

julia> broadcast(+, a, A)
2×3 Array{Float64,2}:
 1.20813  1.82068  1.25387
 1.56851  1.86401  1.67846

julia> b = rand(1,2)
1×2 Array{Float64,2}:
 0.867535  0.00457906

julia> broadcast(+, a, b)
2×2 Array{Float64,2}:
 1.71056  0.847604
 1.73659  0.873631

类似 .+.* 的[点运算符](http://127.0.0.5/@ref man-dot-operators) 等同于 broadcast 调用(除了它们融合两种操作,如下所述)。还有一个 broadcast! 函数来指定一个显式目标(也可以通过 .= 赋值以融合的方式访问它)。事实上,f.(args...) 等价于 broadcast(f, args...),并为广播任何函数提供了方便的语法([dot syntax](http://127.0.0.5/@ref man-vectorized))。嵌套的「点调用」f.(...)(包括调用 .+ 等)[自动融合](http://127.0.0.5/@ref man-dot-operators)到单个 broadcast 调用。

另外,broadcast 并不局限于数组(参见函数文档),对于元组依然有效。对于其他非数组,元组或引用 Ref(除了指针 Ptr) 的参数,视作“标量”。

julia> convert.(Float32, [1, 2])
2-element Array{Float32,1}:
 1.0
 2.0

julia> ceil.((UInt8,), [1.2 3.4; 5.6 6.7])
2×2 Array{UInt8,2}:
 0x02  0x04
 0x06  0x07

julia> string.(1:3, ". ", ["First", "Second", "Third"])
3-element Array{String,1}:
 "1. First"
 "2. Second"
 "3. Third"

实现

Julia 中的基本数组类型是抽象类型 AbstractArray{T,N}。它通过维数 N 和元素类型 T 进行参数化。AbstractVectorAbstractMatrix 是一维和二维情况下的别名。AbstractArray 对象的操作是使用更高级别的运算符和函数定义的,其方式独立于底层存储。这些操作可以正确地被用于任何特定数组实现的回退操作。

AbstractArray 类型包含任何模糊类似的东西,它的实现可能与传统数组完全不同。例如,可以根据请求而不是存储来计算元素。但是,任何具体的 AbstractArray{T,N} 类型通常应该至少实现 size(A)(返回 Int 元组),getindex(A,i) 和 [getindex(A,i1,...,iN)](http://127.0.0.5/@ref getindex);可变数组也应该实现 setindex!。建议这些操作具有几乎为常数的时间复杂性,或严格说来 Õ(1) 复杂性,否则某些数组函数可能出乎意料的慢。具体类型通常还应提供 similar(A,T=eltype(A),dims=size(A)) 方法,用于为 copy 分配类似的数组和其他位于当前数组空间外的操作。无论在内部如何表示 AbstractArray{T,N}T 是由 整数 索引返回的对象类型(A[1, ..., 1],当 A 不为空),N 应该是 size 返回的元组的长度。有关定义自定义 AbstractArray 实现的更多详细信息,请参阅[接口章节中的数组接口导则](http://127.0.0.5/@ref man-interface-array)。

DenseArrayAbstractArray 的抽象子类型,旨在包含元素按列连续存储的所有数组(请参阅[性能建议](http://127.0.0.5/@ref man-performance-tips)中的附加说明)。Array 类型是 DenseArray 的特定实例;VectorMatrix 是在一维和二维情况下的别名。除了所有 AbstractArray 所需的操作之外,很少有专门为 Array 实现的操作;大多数数组库都以通用方式实现,以保证所有自定义数组都具有相似功能。

SubArrayAbstractArray 的特例,它通过与原始数组共享内存而不是复制它来执行索引。 使用view 函数创建 SubArray,它的调用方式与getindex 相同(作用于数组和一系列索引参数)。 view 的结果看起来与 getindex 的结果相同,只是数据保持不变。 view 将输入索引向量存储在 SubArray 对象中,该对象稍后可用于间接索引原始数组。 通过将 @views 宏放在表达式或代码块之前,该表达式中的任何 array [...] 切片将被转换为创建一个 SubArray 视图。

BitArray 是节省空间“压缩”的布尔数组,每个比特(bit)存储一个布尔值。 它们可以类似于 Array{Bool} 数组(每个字节(byte)存储一个布尔值),并且可以分别通过 Array(bitarray)BitArray(array) 相互转换。

「strided」数组存储在内存中,元素以常规偏移量排列,因此元素类型兼容 isbits 的实例可以被传递给期望此内存布局的外部 C 和 Fortran 函数。Strided 数组必须定义一个 strides(A) 方法,该方法为每个维返回一个「strides」元组;提供 stride(A,k) 方法访问该元组中的第 k 个元素。将维数 k 的索引增加 1 应该把 getindex(A,i) 得到的索引 i 增加 stride(A,k)。如果提供了指针转换方法 Base.unsafe_convert(Ptr{T}, A),则内存布局必须以与这些间隔相同的方式对应。DenseArray 是一个非常具体的 strided 数组示例,其中元素是连续排列的,因此它为子类提供了适当的 strides 定义。更多具体的例子可以在 [strided 数组的接口指南](http://127.0.0.5/@ref man-interface-strided-arrays)中找到。StridedVectorStridedMatrix 是许多算是 strided 数组的内置数组类型的便捷别名,允许它们只使用指针和 stride 派发来选择调用专门实现,即使用经高度调试和优化的 BLAS 和 LAPACK 函数。

接下来的例子计算一个大数组中一小部分的 QR 分解,不需要引入任何临时变量。通过正确的维度大小和偏移参数调用合适的 LAPACK 函数。

julia> a = rand(10, 10)
10×10 Array{Float64,2}:
 0.517515  0.0348206  0.749042   0.0979679  …  0.75984     0.950481   0.579513
 0.901092  0.873479   0.134533   0.0697848     0.0586695   0.193254   0.726898
 0.976808  0.0901881  0.208332   0.920358      0.288535    0.705941   0.337137
 0.657127  0.0317896  0.772837   0.534457      0.0966037   0.700694   0.675999
 0.471777  0.144969   0.0718405  0.0827916     0.527233    0.173132   0.694304
 0.160872  0.455168   0.489254   0.827851   …  0.62226     0.0995456  0.946522
 0.291857  0.769492   0.68043    0.629461      0.727558    0.910796   0.834837
 0.775774  0.700731   0.700177   0.0126213     0.00822304  0.327502   0.955181
 0.9715    0.64354    0.848441   0.241474      0.591611    0.792573   0.194357
 0.646596  0.575456   0.0995212  0.038517      0.709233    0.477657   0.0507231

julia> b = view(a, 2:2:8,2:2:4)
4×2 view(::Array{Float64,2}, 2:2:8, 2:2:4) with eltype Float64:
 0.873479   0.0697848
 0.0317896  0.534457
 0.455168   0.827851
 0.700731   0.0126213

julia> (q, r) = qr(b);

julia> q
4×4 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.722358    0.227524  -0.247784    -0.604181
 -0.0262896  -0.575919  -0.804227     0.144377
 -0.376419   -0.75072    0.540177    -0.0541979
 -0.579497    0.230151  -0.00552346   0.781782

julia> r
2×2 Array{Float64,2}:
 -1.20921  -0.383393
  0.0      -0.910506

如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论, 你可以免费 注册 一个本站的账号。
列表为空,暂无数据
    我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站,即表示您同意使用 Cookies 和您的相关数据。
    原文