在Julia中切片和广播多维数组:meshgrid示例

Question

我最近开始通过编写自组织映射的简单实现来学习Julia.我希望用户指定地图的大小和尺寸,这意味着我不能真正使用for循环来处理地图数组,因为我事先并不知道我需要多少层循环.所以我绝对需要广播和切片功能,这些功能适用于任意维度的数组.

现在,我需要构建一个地图索引数组.说我的地图是由大小的数组定义mapsize = (5, 10, 15),我需要构建一个数组indices大小的(3, 5, 10, 15)地方indices[:, a, b, c]应返回[a, b, c].

我来自Python/NumPy背景,其中解决方案已由特定的"函数"给出,mgrid:

indices = numpy.mgrid[:5, :10, :15]
print indices.shape # gives (3, 5, 10, 15)
print indices[:, 1, 2, 3] gives [1, 2, 3]

我没想到朱莉娅在一开始就有这样的功能,所以我转向广播.在NumPy中,广播基于一系列规则,我发现这些规则非常清晰和合乎逻辑.只要每个维度中的大小匹配或其中一个为1,您就可以在同一个表达式中使用不同维度的数组:

(5, 10, 15)   broadcasts to  (5, 10, 15) 
   (10,  1)

(5,  1, 15)   also broadcasts to (5, 10, 15) 
(1, 10,  1)

为此,您还可以使用numpy.newaxis或None轻松地向阵列添加新维度:

array = numpy.zeros((5, 15))
array[:,None,:] has shape (5, 1, 15)

这有助于轻松播放阵列:

A = numpy.arange(5)
B = numpy.arange(10)
C = numpy.arange(15)
bA, bB, bC = numpy.broadcast_arrays(A[:,None,None], B[None,:,None], C[None,None,:])
bA.shape == bB.shape == bC.shape = (5, 10, 15)

使用它,创建indices数组非常简单:

indices = numpy.array(numpy.broadcast_arrays(A[:,None,None], B[None,:,None], C[None,None,:]))
(indices == numpy.mgrid[:5,:10,:15]).all() returns True

一般情况当然有点复杂,但可以使用列表理解和切片来解决:

arrays = [ numpy.arange(i)[tuple([None if m!=n else slice(None) for m in range(len(mapsize))])] for n, i in enumerate(mapsize) ]
indices = numpy.array(numpy.broadcast_arrays(*arrays))

所以回到朱莉娅.我尝试应用相同的理由,最终实现了arrays上面代码列表的等效.由于复合表达式语法,这最终比NumPy对应物更简单:

arrays = [ (idx = ones(Int, length(mapsize)); idx[n] = i;reshape([1:i], tuple(idx...))) for (n,i)=enumerate(mapsize) ]

现在我被困在这里,因为我真的不知道如何将广播应用到我的生成数组列表中...这些broadcast[!]函数要求应用函数f,而我没有.我尝试使用for循环来尝试强制广播:

indices = Array(Int, tuple(unshift!([i for i=mapsize], length(mapsize))...))
for i=1:length(mapsize)
    A[i] = arrays[i]
end

但这给了我一个错误: ERROR: convert has no method matching convert(::Type{Int64}, ::Array{Int64,3})

我这样做是对的吗？我忽略了重要的事情吗？任何帮助表示赞赏.

Answer 1

如果你正在运行julia 0.4,你可以这样做:

julia> function mgrid(mapsize)
           T = typeof(CartesianIndex(mapsize))
           indices = Array(T, mapsize)
           for I in eachindex(indices)
               indices[I] = I
           end
           indices
       end

如果人们可以说,那就更好了

indices = [I for I in CartesianRange(CartesianIndex(mapsize))]

我会调查一下:-).

Answer 2

Julia中的广播几乎都是在NumPy中广播的,所以你应该发现它遵循或多或少相同的简单规则(不确定当不是所有输入具有相同数量的维度时,填充尺寸的方式是否相同但是,因为Julia数组是列专业的.

然而,许多有用的东西,比如newaxis索引,broadcast_arrays还没有实现.(我希望他们愿意.)另请注意,与NumPy相比,Julia中的索引工作方式略有不同:当您在NumPy中留下尾随维度的索引时,其余索引默认为冒号.在朱莉娅,他们可以说是默认为一些.

我不确定你是否真的需要一个meshgrid函数,大多数你想要使用它的东西可以通过使用arrays数组的原始条目和广播操作来完成.meshgrid在matlab中有用的主要原因是因为它在广播中很糟糕.

但是使用该broadcast!函数完成您想要的操作非常简单:

# assume mapsize is a vector with the desired shape, e.g. mapsize = [2,3,4]

N = length(mapsize)
# Your line to create arrays below, with an extra initial dimension on each array
arrays = [ (idx = ones(Int, N+1); idx[n+1] = i;reshape([1:i], tuple(idx...))) for (n,i) in enumerate(mapsize) ]

# Create indices and fill it one coordinate at a time
indices = zeros(Int, tuple(N, mapsize...))
for (i,arr) in enumerate(arrays)
    dest = sub(indices, i, [Colon() for j=1:N]...)
    broadcast!(identity, dest, arr)
end

我必须在条目上添加一个初始单例维度arrays以与轴对齐indices(newaxis在这里很有用......).然后我浏览每个坐标,在相关部分创建一个子阵列(一个视图)indices并填充它.(索引将默认返回Julia 0.4中的子数组,但是现在我们必须sub明确使用).

调用broadcast!仅评估identity(x)=x输入上的标识函数arr=arrays[i],广播到输出的形状.使用该identity功能没有效率损失; broadcast!基于给定函数,参数数量和结果维度数生成专用函数.