Meshgrid function in Julia

Specific Domains Numerics
22

That

is neat! I would have done .* but this trick is much better!

1 Like
23

Is the discussion in the link available in plain-text anywhere? I don’t seem to have access to the page.

24

I just found an easier way to generate the grid data:

x = 1:3;
y = 4:6;
z = 7:9;

xv = getindex.(Iterators.product(x, y, z), 1)  # first.(Iterators.product(x, y, z), 1) is also ok
yv = getindex.(Iterators.product(x, y, z), 2)
zv = getindex.(Iterators.product(x, y, z), 3) # last.(Iterators.product(x, y, z), 3) is also ok

And the result is

julia> getindex.(Iterators.product(x, y, z), 1)
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 1  1  1
 2  2  2
 3  3  3

[:, :, 2] =
 1  1  1
 2  2  2
 3  3  3

[:, :, 3] =
 1  1  1
 2  2  2
 3  3  3

julia> getindex.(Iterators.product(x, y, z), 2)
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 4  5  6
 4  5  6
 4  5  6

[:, :, 2] =
 4  5  6
 4  5  6
 4  5  6

[:, :, 3] =
 4  5  6
 4  5  6
 4  5  6

julia> getindex.(Iterators.product(x, y, z), 3)
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 7  7  7
 7  7  7
 7  7  7

[:, :, 2] =
 8  8  8
 8  8  8
 8  8  8

[:, :, 3] =
 9  9  9
 9  9  9
 9  9  9

which is same as np.mgrid’s result. If you want np.meshgrid style, the first two dimension should be swapped.

4 Likes
25

The easiest way is not to do it at all :wink:

7 Likes
26

To get meshgrid style, x and y shall be swapped as

xv = getindex.(Iterators.product(y, x, z), 2)
yv = getindex.(Iterators.product(y, x, z), 1)

However, if you really need this, by far the easiest and fastest I’ve tested is list comprehension:

function meshgrid(x, y)
   X = [x for _ in y, x in x]
   Y = [y for y in y, _ in x]
   X, Y
end

function meshgrid2(x, y)
   X = getindex.(Iterators.product(y, x), 2)
   Y = getindex.(Iterators.product(y, x), 1)
   X, Y
end

x = 1:1000
y = 1001:2000

julia> @time meshgrid(x,y);
  0.007225 seconds (5 allocations: 15.259 MiB)

julia> @time meshgrid2(x,y);
  0.028644 seconds (9 allocations: 45.777 MiB, 18.41% gc time)

julia> @btime meshgrid($x,$y);
  1.033 ms (4 allocations: 15.26 MiB)

julia> @btime meshgrid2($x,$y);
  4.810 ms (8 allocations: 45.78 MiB)
27

There is now a package that implements such a lazy grid for ngrid (with lazy meshgrid as a one-line variation explained in the docs):

Hopefully this will help a few more people transition from Matlab to Julia :slight_smile:

9 Likes
28

function meshgrid(n,L) #n= number of grids, L= size of interval
a=transpose(range(0,L,length=n))
b=repeat(a,n,1)
return b
end

29

EDIT: I had a comparison with and without using Meshgrid, but I removed it because there was a mistake in the code and the results were erroneous (on top of using time in the global scope)

The part about LazyGrids being faster seemed to be correct though.

EDIT: using LazyGrids (https://github.com/JuliaArrays/LazyGrids.jl)
made it even faster
0.000021 seconds (15 allocations: 928 bytes)

30

Don’t time in the global scope. You’re measuring compilation time.

31

change your name from jlchan to jichan

1 Like
32

You mean Oji-chan?

2 Likes
33

:laughing: yup

1 Like
34

I was just wondering the same thing and I used something like

julia> x_i = 0:1:3
0:1:3

julia> y_i = 0:1:3
0:1:3

julia> x_co = [x for x in x_i, y in y_i]
4×4 Matrix{Int64}:
 0  0  0  0
 1  1  1  1
 2  2  2  2
 3  3  3  3

julia> y_co = [y for x in x_i, y in y_i]
4×4 Matrix{Int64}:
 0  1  2  3
 0  1  2  3
 0  1  2  3
 0  1  2  3

I guess it’s not needed but still just added if someone needs.

35

You’re correct that it isn’t needed and actually it isn’t a good idea when there are better alternatives in Julia. See this answer its follow-up responses for more detailed discussion.

1 Like