Utilizing ndgrid/meshgrid functionality in Julia

Question

I'm trying to find functionality in Julia similar to MATLAB's meshgrid or ndgrid. I know Julia has defined ndgrid in the examples but when I try to use it I get the following error.

UndefVarError: ndgrid not defined

Anyone know either how to get the builtin ndgrid function to work or possibly another function I haven't found or library that provides these methods (the builtin function would be preferred)? I'd rather not write my own in this case.

Thanks!

Pronunciamento

The question is actually incomplete - much like asking "how would you translate the English word 'of' into German". The case is that in most code in matlab where ndgrid is ideomatic, the Julia code would rely on dot-broadcasting. Try expanding the question to describe what you want to achieve. – Assess

VectorizedRoutines.jl has meshgrid and ndgrid functions. However, I think that the proper way to do this would be by the creation of a lazy operator as output. Never got to it though. – Empty

How are the function definitions in the linked file run? using include("$JULIA_HOME/../examples/ndgrid.jl") or some such statement? – Enchantment

I've updated the wording to hopefully clarify my question. I think @ChrisRackauckas may have answered it for now but if there was a way to use the builtin function and not an outside lib that would be preferred. Thanks so much for the help so far! – Pronunciamento

Why built in an not a package? What's the difference in Julia? – Empty

Honestly I just thought it would be easier to distribute my code. Not a huge deal. I'm new to Julia so I hoping someone would be like "oh yeah you just left this out...". – Pronunciamento

@ChrisRackauckas I ended up using your VectorizedRountines.jl package. Works great. Feel free to add that as the answer to the question. Thanks! – Pronunciamento

Answer 1

We prefer to avoid these functions, since they allocate arrays that usually aren't necessary. The values in these arrays have such a regular structure that they don't need to be stored; they can just be computed during iteration. For example, one alternative approach is to write an array comprehension:

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

Or, you can write for loops, or iterate over a product iterator:

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

Ensoll

Wouldn't it be nice to have an AbstractArray type that computes the values as necessary without allocations just like linspace though? Essentially an n-dimensional generalization of linspace. – Encarnalize

Answer 2

I do find sometimes it's convenient to use some function like meshgrid in numpy. It's easy to do it with list comprehension:

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

e.g.

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

now

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

However, I did not find this makes the code run faster than using iteration. Here's what I mean:

After defining

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

I did benchmark on the following two functions

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

Here are the benchmark results:

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

The meshgrid method is both significantly slower and taking more memory. Any Julia expert knows why? I understand Julia is a compiling language unlike Python so iterations won't be slower than vectorization, but I don't understand why vector(array) calculation is many times slower than iteration. (For bigger N this difference is even larger.)

Edit: After reading this post, I have the following updated version of the 'meshgrid' method. The idea is to not create a meshgrid beforehand, but to do it in the calculation via Julia's powerful elementwise array operation:

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

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

The trick here is the .+ between a size-M column array and a size-N row array which returns a M-by-N array. It does the 'meshgrid' for you. This function is nearly 3 times faster then fun1, albeit not as fast as fun2.

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

Gower

Two things. First, fun1 should be mean(sqrt.(X.*X .+ Y.*Y)) as this saves an array copy. That said, the bigger reason is simply that making arrays is more expensive than not making arrays. Also, in this case your arrays are non-const globals, which is often very slow. – Jounce

After removing the 'return' and declaring X and Y as const, the computation time for func1 is still 8.114 ms, only 0.2 ms faster. So, array computations are simply slow in Julia. My takeaway from this is to code Julia more like C and less like Python or MATLAB, i.e. use iterations instead of arrays. – Gower

Answer 3

Above, @ChrisRackauckas suggests that the "proper way" to do this is with a lazy operator but he hadn't gotten around to it.

There is now a registered packaged with lazy ndgrid in it: https://github.com/JuliaArrays/LazyGrids.jl

It is more general than the version in VectorizedRoutines.jl because it can handle vectors with different types, e.g., ndgrid(1:3, Float16[0:2], ["x", "y", "z"]).

There are Literate.jl examples in the docs that show the lazy performance is pretty good.

Of course lazy meshgrid is just one step away:

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

Lengthwise

Answer 4

Given two axes

x = 1:4
y = 1:3

one can generated a meshgrid with

X = x' .+ 0y
Y = 0*x' .+ y

[Multiplying x with 0 requires the asterics, since 0x is attempted to parse as a hex number (e.g. 0x0A) and would yield invalid numeric constant]
This code yields

3×4 Matrix{Int64}:
 1  2  3  4
 1  2  3  4
 1  2  3  4
3×4 Matrix{Int64}:
 1  1  1  1
 2  2  2  2
 3  3  3  3

To generate a flattened meshgrid, one can use

X = repeat(x, length(y))
Y = repeat(y, inner=length(x))

# [1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4]
# [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3]

Or pairwise

julia> [(xi, yi) for yi=y for xi=x]
12-element Vector{Tuple{Int64, Int64}}:
 (1, 1)
 (2, 1)
 (3, 1)
 (4, 1)
 (1, 2)
 (2, 2)
 (3, 2)
 (4, 2)
 (1, 3)
 (2, 3)
 (3, 3)
 (4, 3)

Sb