I am also confused. I use this for contour
in PyPlot, where the help files tell me to use meshgrid. How can I avoid meshgrid without finding a different plotting package?
You don’t need a meshgrid for contour plotting with PyPlot — you can pass 1d arrays for the axes and it knows how to broadcast them, and you write 2d functions of the axes variables by broadcast operations as well:
using PyPlot
x = range(0,2,length=100)' # note ': this is a row vector
y = range(0,4,length=200)
z = @. sin(x) * cos(y) # broadcasts to 2d array
contour(x, y, z)
Thanks. I’m even older than @jlchan and transcribed my stuff from (no joke) 20 year old matlab.
This is the corresponding Matlab code, btw:
x = linspace(0, 2, 100);
y = linspace(0, 4, 200).';
z = sin(x) .* cos(y); % broadcasts to 2d array
contour(x, y, z)
I have seen the error in my ways and will go forth and sin in a different way in the future.
Remember to cos as well!
I will do that as I sit in the sun and work on my tan.
Instead of sin
ning, sec
and ye shall find!
That
is neat! I would have done .* but this trick is much better!
Is the discussion in the link available in plain-text anywhere? I don’t seem to have access to the page.
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.
The easiest way is not to do it at all
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)
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
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
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)
Don’t time in the global scope. You’re measuring compilation time.
change your name from jlchan to jichan
You mean Oji-chan?
yup