注意：更新于 2018 年 12 月，适用于 Julia 1.1

注意：更新于 2020 年 4 月，为清晰起见

Julia 使编写优雅高效的多维算法变得容易。新的功能建立在两个基础之上：一个名为 CartesianIndices 的迭代器，以及复杂的数组索引机制。在解释之前，我想强调一下，开发这些功能是一个协作过程，大部分工作是由 Matt Bauman (@mbauman)、Jutho Haegeman (@Jutho) 和我 (@timholy) 完成的。

这些迭代器看似简单：仅仅几个原则就能在编写多维算法时带来强大的力量。然而，就像许多简单的概念一样，其含义可能需要一段时间才能理解。也有可能将这些技术与 Base.Cartesian 混淆，而 Base.Cartesian 是一种完全不同的（更痛苦的）解决同一问题的方案。仍然有一些情况需要使用 Base.Cartesian，但对于许多问题来说，这些新功能代表着一种简化得多的方法。

让我们用一个从 手册 中摘取的示例的扩展来介绍这些迭代器。

eachindex、CartesianIndex 和 CartesianIndices

有两种推荐的迭代 AbstractArray 中元素的“默认”方法：如果您不需要与每个元素相关联的索引，则可以使用

for a in A    # A is an AbstractArray
    # Code that does something with the element a
end

如果您还需要索引，则使用

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

在某些情况下，此循环的第一行扩展为 for i = 1:length(A)，而 i 只是一个整数。但是，在其他情况下，它将扩展为等效于

for i in CartesianIndices(A)
    # i is now a CartesianIndex
    # Code that does something with i and/or A[i]
end

您可以使用以下代码自行查看它的作用

julia> A = rand(3,2)

julia> for i in CartesianIndices(A)
          @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)

CartesianIndex{N} 表示一个 N 维索引。CartesianIndex 基于元组，实际上您可以使用 Tuple(i) 访问底层元组。

CartesianIndices 就像一个 CartesianIndex 值的数组

julia> iter = CartesianIndices(A)
3×2 CartesianIndices{2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}}:
 CartesianIndex(1, 1)  CartesianIndex(1, 2)
 CartesianIndex(2, 1)  CartesianIndex(2, 2)
 CartesianIndex(3, 1)  CartesianIndex(3, 2)

julia> supertype(typeof(iter))
AbstractArray{CartesianIndex{2},2}

因此，iter[2,2] 和 iter[5] 都返回 CartesianIndex(2, 2)；实际上，后者是将 线性索引 转换为多维笛卡尔索引的推荐方法。

但是，在内部，iter 只是一个围绕每个维度的 axes 范围的包装器

julia> iter.indices
(Base.OneTo(3), Base.OneTo(2))

因此，在许多应用程序中，这些对象的创建和使用几乎没有或根本没有开销。

您可以手动构造它们：例如，

julia> CartesianIndices((-7:7, 0:15))
15×16 CartesianIndices{2,Tuple{UnitRange{Int64},UnitRange{Int64}}}:
# remaining output suppressed

对应一个将沿着第一个维度遍历 -7:7、沿着第二个维度遍历 0:15 的迭代器。

eachindex 比 for i = 1:length(A) 更推荐的原因之一是，某些 AbstractArray 无法使用线性索引有效地进行索引；相反，更广泛的对象类别可以有效地使用多维迭代器进行索引。（一般来说，SubArrays 是一个主要示例。）eachindex 被设计为为给定数组类型选择最有效的迭代器。您甚至可以使用

for i in eachindex(A, B)
    ...

来增加 i 有效地访问 A 和 B 的可能性。（使用 eachindex 的第二个原因是一些数组的索引不从 1 开始，但这是另一篇 博客文章 的主题。）

正如我们将在下面看到的那样，这些迭代器还有另一个目的：无论底层数组是否具有有效的线性索引，多维迭代在编写算法时都可以成为强大的盟友。这篇博客文章的其余部分将重点介绍后一种应用。

使用 CartesianIndex 迭代器编写多维算法

多维箱形滤波器

假设我们有一个多维数组 A，并且我们希望计算每个元素周围 3x3x... 块的 “移动平均值”。从任何给定的索引位置，我们将希望沿着每个维度对偏移 -1:1 的区域进行求和。当然，边缘位置必须进行特殊处理，以避免超出数组的边界。

在许多语言中，编写此概念上简单的算法的通用（N 维）实现会比较痛苦，但在 Julia 中这非常容易

function boxcar3(A::AbstractArray)
    out = similar(A)
    R = CartesianIndices(A)
    Ifirst, Ilast = first(R), last(R)
    I1 = oneunit(Ifirst)
    for I in R
        n, s = 0, zero(eltype(out))
        for J in max(Ifirst, I-I1):min(Ilast, I+I1)
            s += A[J]
            n += 1
        end
        out[I] = s/n
    end
    out
end

（请注意，此示例仅适用于 Julia 版本 1.1 及更高版本。）

让我们逐行分析

• out = similar(A) 为输出分配内存。在“真实”实现中，您需要对输出的元素类型更加小心（如果输入数组的元素类型是 Int 会怎么样？），但为了简单起见，我们在这里简化了一些内容。

• R = CartesianIndices(A) 为数组创建迭代器。假设 A 的索引从 1 开始，则范围从 CartesianIndex(1, 1, 1, ...) 到 CartesianIndex(size(A,1), size(A,2), size(A,3), ...)。我们不使用 eachindex，因为我们无法确定它是否会返回一个 CartesianIndices 迭代器，而在这里我们明确需要一个。

• Ifirst = first(R) 和 Ilast = last(R) 分别返回迭代范围的下界（CartesianIndex(1, 1, 1, ...)）和上界（CartesianIndex(size(A,1), size(A,2), size(A,3), ...)）。我们将使用它们来确保我们永远不会访问 A 的越界元素。

• I1 = oneunit(Ifirst) 创建一个与 Ifirst 维度相同的全 1 CartesianIndex。我们将在算术运算中使用它来定义感兴趣区域。

• for I in R：在这里，我们遍历 R 的每个条目，对应于 A 和 out。

• n = 0 和 s = zero(eltype(out)) 初始化累加器。s 将保存相邻值的总和。n 将保存使用的相邻元素的数量；在大多数情况下，在循环结束后，我们将有 n == 3^N，但对于边缘点，有效相邻元素的数量会更小。

• for J in max(Ifirst, I-I1):min(Ilast, I+I1) 可能是算法中最“巧妙”的一行。I-I1 是一个 CartesianIndex，它沿着每个维度都比 I 低 1，而 I+I1 则高 1。但是，当 I 表示一个边缘点时，I-I1 或 I+I1（或两者）都可能越界。max(Ifirst, I-I1) 确保 J 的每个坐标都大于或等于 1，而 min(Ilast, I+I1) 确保 J[d] <= size(A,d)。

  将这两个元素与冒号组合起来，Ilower:Iupper，创建一个用作迭代器的 CartesianIndices 对象。

• 内部循环累积 s 中的总和以及 n 中访问的相邻元素的数量。

• 最后，我们将平均值存储在 out[I] 中。

此实现不仅简单，而且非常健壮：对于边缘点，它计算其所有可用最近邻居的平均值。即使对于某些维度 d，size(A, d) < 3 也能正常工作；我们不需要对 A 的大小进行任何错误检查。

计算约简

作为第二个示例，考虑多维约简的实现。约简接受一个输入数组，并返回一个尺寸更小的数组（或标量）。一个经典的示例是沿着数组的特定维度进行求和：对于三维数组，您可能希望沿着维度 2 进行求和，保留维度 1 和 3 不变。

核心算法

编写此算法的一种有效方法要求输出数组 B 由调用者预先分配（稍后我们将看到如何以编程方式分配 B）。例如，如果输入 A 的大小为 (l,m,n)，那么在仅沿着维度 2 进行求和时，输出 B 的大小将为 (l,1,n)。

在给定此设置的情况下，实现非常简单

function sumalongdims!(B, A)
    # It's assumed that B has size 1 along any dimension that we're summing,
    # and otherwise matches A
    fill!(B, 0)
    Bmax = last(CartesianIndices(B))
    for I in CartesianIndices(A)
        B[min(Bmax,I)] += A[I]
    end
    B
end

此算法背后的关键思想封装在单个语句 B[min(Bmax,I)] 中。对于我们三维示例，其中 A 的大小为 (l,m,n)，B 的大小为 (l,1,n)，内部循环本质上等效于

B[i,1,k] += A[i,j,k]

因为 min(1,j) = 1。

包装器，以及使用函数屏障处理类型不稳定性

作为用户，您可能更喜欢类似于 sumalongdims(A, dims) 的接口，其中 dims 指定您要沿着其求和的维度。dims 可以是单个整数，例如我们上面示例中的 2，也可以是（如果您希望一次沿着多个维度求和）元组或 Vector{Int}。这实际上是 sum(A; dims=dims) 中使用的接口；在这里，我们希望编写我们自己的（稍微简单一些的）实现。

包装器的一种可能的简化实现如下所示

function sumalongdims(A, dims)
    sz = [size(A)...]
    sz[[dims...]] .= 1
    B = Array{eltype(A)}(undef, sz...)
    sumalongdims!(B, A)
end

显然，此简单的实现跳过了所有相关的错误检查。但是，这里我想要探讨的主要点是，B 的分配实际上是 不可推断的：sz 是一个 Vector{Int}，特定 Vector{Int} 的长度（元素数量）没有被类型本身编码，因此 B 的维度无法推断。

现在，我们可以通过多种方式解决这个问题，例如通过对结果进行注释

B = Array{eltype(A)}(undef, sz...)::Array{eltype(A),ndims(A)}

或通过使用可推断的实现

function sumalongdims(A, dims)
    sz = ntuple(i->i ∈ dims ? 1 : size(A, i), Val(ndims(A)))
    B = Array{eltype(A)}(undef, sz...)
    sumalongdims!(B, A)
end

但是，在这里我们想强调的是，这种设计——从 sumalongdims 中分离出 sumalongdims!——通常可以缓解推断问题最严重的影响。这种技巧，使用 函数调用将性能关键步骤与潜在的类型不稳定的前置步骤隔离开来，有时被称为引入函数屏障。它允许 Julia 的编译器生成 sumalongdims! 的良好优化版本，即使 B 的中间类型未知。

一般来说，在编写多维代码时，您应该确保主迭代与类型不稳定的前置步骤处于不同的函数中。（在早期版本的 Julia 中，您可能会看到使用 @noinline 注释的内核函数，以防止内联器将两者重新组合在一起，但对于较新的版本的 Julia，这不再必要。）

当然，在这个示例中，将它制作成一个独立的函数还有第二个动机：如果此计算是您将要多次重复的计算，那么重复使用同一个输出数组可以减少代码中的内存分配量。

沿着指定维度进行滤波（利用多个索引）

最后一个示例说明了一个重要的新要点：当您索引数组时，可以随意混合 CartesianIndex 和整数。为了说明这一点，我们将编写一个 指数平滑滤波器。实现此类滤波器的有效方法是使平滑后的输出值 s[i] 依赖于当前输入 x[i] 和之前的滤波值 s[i-1] 的组合；在一维中，您可以将其写为

function expfilt1!(s, x, α)
    0 < α <= 1 || error("α must be between 0 and 1")
    s[1] = x[1]
    for i = 2:length(x)
        s[i] = α*x[i] + (1-α)*s[i-1]
    end
    s
end

这将导致时间尺度为 1/α 的近似指数衰减。

这里，我们希望实现该算法，以便它可以用来沿任何选定维度对数组进行指数滤波。再次强调，该实现出奇地简单。

function expfiltdim(x, dim::Integer, α)
    s = similar(x)
    Rpre = CartesianIndices(size(x)[1:dim-1])
    Rpost = CartesianIndices(size(x)[dim+1:end])
    _expfilt!(s, x, α, Rpre, size(x, dim), Rpost)
end

function _expfilt!(s, x, α, Rpre, n, Rpost)
    for Ipost in Rpost
        # Initialize the first value along the filtered dimension
        for Ipre in Rpre
            s[Ipre, 1, Ipost] = x[Ipre, 1, Ipost]
        end
        # Handle all other entries
        for i = 2:n
            for Ipre in Rpre
                s[Ipre, i, Ipost] = α*x[Ipre, i, Ipost] + (1-α)*s[Ipre, i-1, Ipost]
            end
        end
    end
    s
end

再次注意函数屏障技术的应用。在核心算法（_expfilt!）中，我们的策略是使用两个CartesianIndex迭代器，Ipre 和 Ipost，其中第一个涵盖维度1:dim-1，第二个涵盖dim+1:ndims(x)；过滤维度dim由一个整数索引i单独处理。由于过滤维度由一个整数输入指定，因此无法推断每个索引元组Ipre 和 Ipost 中有多少个条目。因此，我们在算法的类型不稳定部分计算CartesianIndices，然后将其作为参数传递给核心例程_expfilt!。

使这种实现成为可能的是，我们可以将x 索引为x[Ipre, i, Ipost]。请注意，提供的索引总数为(dim-1) + 1 + (ndims(x)-dim)，这正好是ndims(x)。通常，在对 Julia 中的AbstractArray进行索引时，可以提供任何整数和CartesianIndex索引的组合。

AxisAlgorithms 包广泛地使用了这些技巧，并反过来为像 Interpolations 这样的高性能包提供了核心支持，这些包需要多维计算。

其他问题

缓存效率

值得注意的是，到目前为止还没有提到的一点：这里的所有示例都是相对缓存效率的。这是编写高效代码时要注意的关键属性。特别是，julia 数组按从第一个到最后一个维度的顺序存储（对于矩阵，是“列优先”顺序），因此应该从最后一个到第一个维度嵌套迭代。例如，在上面的过滤示例中，我们小心地按照以下顺序进行迭代

for Ipost ...
    for i ...
        for Ipre ...
            x[Ipre, i, Ipost] ...

以便按内存顺序遍历x。

广播

CartesianIndexes不是可广播的

julia> I = CartesianIndex(2, 7)
CartesianIndex(2, 7)

julia> I .+ 1
ERROR: iteration is deliberately unsupported for CartesianIndex. Use `I` rather than `I...`, or use `Tuple(I)...`
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] iterate(::CartesianIndex{2}) at ./multidimensional.jl:154
...

当要执行广播运算时，只需提取底层的元组

julia> Tuple(I) .+ 1
(3, 8)

如果需要，可以将其打包回CartesianIndex，或者直接使用它（带展开）进行索引。编译器优化了所有这些操作，因此这种方式创建对象没有实际“成本”。

为什么不允许迭代？原因之一是为了支持以下情况

julia> R = CartesianIndices((1:3, 1:3))
3×3 CartesianIndices{2,Tuple{UnitRange{Int64},UnitRange{Int64}}}:
 CartesianIndex(1, 1)  CartesianIndex(1, 2)  CartesianIndex(1, 3)
 CartesianIndex(2, 1)  CartesianIndex(2, 2)  CartesianIndex(2, 3)
 CartesianIndex(3, 1)  CartesianIndex(3, 2)  CartesianIndex(3, 3)

julia> R .+ CartesianIndex(2, 17)
3×3 CartesianIndices{2,Tuple{UnitRange{Int64},UnitRange{Int64}}}:
 CartesianIndex(3, 18)  CartesianIndex(3, 19)  CartesianIndex(3, 20)
 CartesianIndex(4, 18)  CartesianIndex(4, 19)  CartesianIndex(4, 20)
 CartesianIndex(5, 18)  CartesianIndex(5, 19)  CartesianIndex(5, 20)

基本思想是CartesianIndex(2, 17) 在任何地方都必须像一对标量索引一样工作；因此，CartesianIndex 必须被视为一个单独的（标量）实体，而不是一个容器本身。

总结

如现在希望清楚的那样，Julia 的优雅迭代器消除了编写通用多维算法的大部分痛苦。这里的示例只是触及了表面，但底层原理非常简单；希望这些示例能让你更容易编写自己的算法。