在 Julia 中使用 ndgrid/meshgrid 功能答案

【问题标题】：Utilizing ndgrid/meshgrid functionality in Julia在 Julia 中使用 ndgrid/meshgrid 功能
【发布时间】：2021-07-12 10:35:27
【问题描述】：

我试图在 Julia 中找到类似于 MATLAB 的 meshgrid 或 ndgrid 的功能。我知道 Julia 已经在 the examples 中定义了 ndgrid，但是当我尝试使用它时，出现以下错误。

UndefVarError: ndgrid 未定义

任何人都知道如何让内置ndgrid 函数工作，或者可能是我还没有找到的另一个函数或提供这些方法的库（首选内置函数）？在这种情况下，我宁愿不自己写。

谢谢！

【问题讨论】：

这个问题实际上是不完整的——就像问“你如何将英语单词'of'翻译成德语”。情况是，在 ndgrid 是 ideomatic 的 matlab 中的大多数代码中，Julia 代码将依赖于点广播。尝试扩展问题以描述您想要实现的目标。
VectorizedRoutines.jl 具有 meshgrid 和 ndgrid 函数。但是，我认为正确的方法是创建一个惰性运算符作为输出。不过从来没有遇到过。
链接文件中的函数定义是如何运行的？使用include("$JULIA_HOME/../examples/ndgrid.jl") 或类似的声明？
我已经更新了措辞，希望能澄清我的问题。我认为@ChrisRackauckas 现在可能已经回答了这个问题，但是如果有办法使用内置函数而不是首选的外部库。非常感谢到目前为止的帮助！
为什么要内置一个非包？ Julia 有什么不同？

标签： julia

【解决方案1】：

我们宁愿避免使用这些函数，因为它们分配的数组通常是不必要的。这些数组中的值具有如此规则的结构，因此不需要存储它们；它们只能在迭代期间计算。例如，另一种方法是编写数组推导：

julia> [ 10i + j for i=1:5, j=1:5 ]
5×5 Array{Int64,2}:
 11  12  13  14  15
 21  22  23  24  25
 31  32  33  34  35
 41  42  43  44  45
 51  52  53  54  55

或者，您可以编写 for 循环，或迭代 product 迭代器：

julia> collect(Iterators.product(1:2, 3:4))
2×2 Array{Tuple{Int64,Int64},2}:
 (1, 3)  (1, 4)
 (2, 3)  (2, 4)

【讨论】：

如果有一个 AbstractArray 类型可以根据需要计算值而不像 linspace 那样进行分配，这不是很好吗？本质上是 linspace 的 n 维泛化。

【解决方案2】：

我确实发现有时在 numpy 中使用像 meshgrid 这样的函数很方便。使用列表理解很容易做到这一点：

function meshgrid(x, y)
    X = [i for i in x, j in 1:length(y)]
    Y = [j for i in 1:length(x), j in y]
    return X, Y
end

例如

x = 1:4
y = 1:3
X, Y = meshgrid(x, y)

现在

julia> X
4×3 Array{Int64,2}:
 1  1  1
 2  2  2
 3  3  3
 4  4  4
julia> Y
4×3 Array{Int64,2}:
 1  2  3
 1  2  3
 1  2  3
 1  2  3

但是，我没有发现这会使代码运行得比使用迭代更快。这就是我的意思：

定义后

x = 1:1000
y = x
X, Y = meshgrid(x, y)

我对以下两个函数做了基准测试

using Statistics

function fun1()
    return mean(sqrt.(X.*X + Y.*Y))
end

function fun2()
    sum = 0.0
    for i in 1:1000
        for j in 1:1000
            sum += sqrt(i*i + j*j)
        end
    end
    return sum / (1000*1000)
end

以下是基准测试结果：

julia> @btime fun1()
  8.310 ms (19 allocations: 30.52 MiB)
julia> @btime run2()
  1.671 ms (0 allocations: 0 bytes)

meshgrid 方法速度慢得多，而且占用更多内存。 任何 Julia 专家都知道为什么吗？我了解 Julia 是一种与 Python 不同的编译语言，因此迭代不会比向量化慢，但我不明白为什么向量（数组）计算比迭代慢很多倍。（对于更大的 N，这个差异甚至更大。）

编辑： 阅读this post 后，我有以下更新版本的“meshgrid”方法。这个想法不是事先创建网格网格，而是通过 Julia 强大的元素数组操作在计算中完成：

x = collect(1:1000)
y = x'

function fun1v2()
    mean(sqrt.(x .* x .+ y .* y))
end

这里的技巧是大小为 M 的列数组和大小为 N 的行数组之间的 .+ 返回一个 M×N 数组。它为你做“网格”。这个函数比fun1快了近3倍，虽然没有fun2快。

julia> @btime fun1v2()
  3.189 ms (24 allocations: 7.63 MiB)
765.8435104896155

【讨论】：

两件事。首先，fun1 应该是 mean(sqrt.(X.*X .+ Y.*Y))，因为这样可以保存数组副本。也就是说，更大的原因仅仅是制造阵列比不制造阵列更昂贵。此外，在这种情况下，您的数组不是const 全局变量，这通常非常慢。
去掉'return'并将X和Y声明为const后，func1的计算时间仍然是8.114 ms，只快了0.2 ms。因此，Julia 中的数组计算很慢。我的收获是将 Julia 编码得更像 C 而不是 Python 或 MATLAB，即使用迭代而不是数组。

【解决方案3】：

在上面，@ChrisRackauckas 建议“正确的方法”是使用懒惰的操作员，但他还没有解决这个问题。

现在有一个注册包，里面有惰性ndgrid： https://github.com/JuliaArrays/LazyGrids.jl

它比中的版本更通用 VectorizedRoutines.jl 因为它可以处理不同类型的向量，例如， ndgrid(1:3, Float16[0:2], ["x", "y", "z"]).

文档中有Literate.jl 示例表明惰性性能非常好。

当然懒meshgrid就差一步了：

meshgrid(y,x) = (ndgrid_lazy(x,y)[[2,1]]...,)

【讨论】：