Julia: add Matlab-like meshgrid function

[stevengj]

I've always found Matlib's meshgrid function to be very useful for plotting and for working with data on Cartesian grids.

Would be nice to have this in Julia, and is trivial to implement. It could even return an AbstractArray type that avoided allocating the full matrices of (redundant) numbers, although to be honest it's probably not worth it (if you look at cases where meshgrid is used, in my experience the memory usage is rarely an issue — you always do some computations with the meshgrid that requires comparable storage, so a compressed format for the meshgrid will at best save you a small constant factor in memory usage).

[ViralBShah]

I believe we have ndgrid. There is also quite a bit of discussion on this in older issues.

[stevengj]

It looks like ndgrid is just in examples, which is basically useless (most users will never see it).

NumPy also calls it meshgrid. (And also has an mgrid function.) I never realized that Matlab had an ndgrid too, ugh!

[stevengj]

I don't understand the opposition to meshgrid (or ndgrid) in #3026. What is the concise alternative to, for example:

using PyPlot
x, y = meshgrid(linspace(0,2pi,50), linspace(0,2pi,50))
contour(x, y, sin(x + cos(x.*y)))

@JeffBezanson? @timholy?

[timholy]

Mostly I'm opposed to having both meshgrid and ndgrid. Having two functions that do the same thing except for transposition in the first two coordinates is silly and confusing. To me it's easier to remember an API like ndgrid's, which takes inputs of length m, n, p, ... and generates an output of size m-by-n-by-p-by..., than it is meshgrid's (which generates an output of size n-by-m-by-p-by...).

Even more radically, in a (non-realistic) world of pure Julia one wonders if either one is needed. Given that we have broadcasting and multiple dispatch I don't see much need to have a function that takes two 1d objects and blows them up into a 2d object. I.e., in Julia your example should be

y = linspace(0,2pi,50)
x = y'
contour(x, y, sin(x .+ cos(x .* y)))

However, if one is interfacing with outside packages then one needs such functions. And even within Julia users may want both. I just think it's ugly.

[stevengj]

Hmm, something like that _should_ work with matplotlib, now that I look at their contour docs. But I get:

no method ##*_inner!#201((Int64,Int64),(Array{Float64,2},Array{Float64,1}),Array{Int64,1},Array{Union(Uint64,Int64),2},Array{Float64,2})
at In[83]:3
 in ##broadcast_T_*#199 at broadcast.jl:178
 in .* at broadcast.jl:265

is something wrong with the broadcasting functions?

[timholy]

I've seen that once myself, but when I tried following it up on a later occasion I couldn't reproduce it. It must depend on particulars in a way that I didn't immediately discover.

[StefanKarpinski]

@stevengj, if you have a reproducible case, please do open an issue for it recording how you produced the issue.

[JeffBezanson]

There are implementations of ndgrid and meshgrid in examples/ndgrid.jl. I agree ndgrid makes much more sense. But I've always disliked the way they store O(d) information in O(n^d) space.

[stevengj]

As I mentioned above, the space wastage is in practice only a constant factor, since you always use these arrays to compute something else.

[JeffBezanson]

True, but reading the parameters from memory instead of computing them with a couple instructions still can't be good. I don't have numbers at hand to back that up, but there is a palpable silliness to it.

[stevengj]

You are still arguing about constant factors. As soon as one is using a sequence of vector operations rather than folding everything into a single loop, one has already lost the constant-factor battle. Things like meshgrid are for convenience, not for optimizing inner loops.

[JeffBezanson]

I suppose there is no problem with copying the ~30 lines of code from examples to Base.

The matlab API is designed to lose the constant-factor battle. If functions like contour accepted range objects and a function, there would be a fighting chance of folding everything into a single loop.

[JeffBezanson]

Would it be ok to add just ndgrid?

[stevengj]

I suppose, although the transposition is useful for plotting.

[StefanKarpinski]

Perhaps a gridplot function would be in order that takes things in the expected order?

[JeffBezanson]

If we were allowed to change the plotting api, this wouldn't come up. This
is for calling matplotlib.
On Aug 30, 2013 11:28 AM, "Stefan Karpinski" [email protected]
wrote:

Perhaps a gridplot function would be in order that takes things in the
expected order?

—
Reply to this email directly or view it on GitHubhttps://github.com/JuliaLang/julia/issues/4093#issuecomment-23569057
.

[stevengj]

I can change the plotting API if needed, if we are only talking about matplotlib, if you are calling it through PyPlot. So let's just stick with ndgrid for now and revisit the issue later if needed.

[stevengj]

(I've updated PyPlot to consistently accept 1d arrays or row/column vectors for the first two arguments in contour, surface, etc. plots. It turned out that Matplotlib was already doing this for some functions, but it wasn't consistent.)

[tknopp]

One use case for meshgrid is drawing shapes in images. Consider for instance drawing a circle into an array by:

x = linspace(-1,1,256)
xx,yy = meshgrid( x, x )
im = xx.^2 + yy.^2 < (0.2)^2

What would be the Julian way to do this? As @stevengj pointed out, if one worries about the performance / memory consumption, one could write a dedicated type that calculates the values on the fly by implementing getindex. But as an initial implementiation, the "memory hungry" meshgrid function seems to me just fine.

[stevengj]

Just use broadcasting operations:

x = linspace(-1,1, 200)
y = x'
im = x.^2 .+ y.^2 .< 0.2^2

(My need for the meshgrid feature has decreased a lot because I just made sure that the PyPlot wrappers support broadcasting of 1d row/column arrays.)

[tknopp]

A thanks! Might be something for the Matlab FAQ in the manual

[ihnorton]

(please re-open and tag decision or up-for-grabs if this needs further discussion...)

[pluskid]

Sometimes I still need to generate grid points explicitly, because I need to manipulate those points via some procedures that does not support this kind of broadcasting. It would be great if a convenient meshgrid like function is available.

[ViralBShah]

Won't a comprehension do it? Slightly longer at best. Maybe you could have it as a convenience function in your own code?

[pluskid]

@ViralBShah Can you show me how to do it with comprehension? I'm currently doing it with a simple and dumb way like this:

  a = linspace(0, pi, n)
  b = linspace(0, 2pi, 2n)
  grid_a = vec(broadcast((x,y) -> x, a, b'))
  grid_b = vec(broadcast((x,y) -> y, a, b'))

[IainNZ]

grid_a = [i for i in a, j in b]
grid_b = [j for i in a, j in b]

[pluskid]

@IainNZ Thanks!

[lkuper]

Any chance of releasing the examples/ndgrid.jl implementation of ndgrid as part of a library? (I'm asking because we have a ParallelAccelerator example that relies on it, and to release that example it seems like we'd need to copy that file into our own project and make a note in our license file that we are redistributing a file that's part of Julia.)

[stevengj]

Wouldn't it be better to rewrite your ParallelAccelerator example to use broadcasting?

[lkuper]

@stevengj It might be! We're benchmarking it against a Matlab version that calls Matlab's meshgrid (I meant to say meshgrid above, not ndgrid), and I'm trying to strike a balance between doing what is idiomatic for Julia and doing something that would be familiar to a Matlab programmer.

[SimonEnsemble]

Yes, not having a meshgrid function has caused difficulty multiple times for me (a user) when trying to make visualizations.

[pengwyn]

I also have come across a different use case. There appears to be no broadcast equivalent to pmap, so it is necessary to grid my vectors across into same-sized n-dimensional arrays. It would be much nicer to instead write pmap(func, x', y).

Alternatively, a different function could solve my problem. Maybe #10438 is relevant here?

[johnnychen94]

We could implement a helper like np.ix_

function xvec(vs::AbstractVector...)
    n = length(vs)
    map(enumerate(vs)) do (i, v)
        sz = ntuple(x->x==i ? length(v) : 1, n)
        reshape(v, sz)
    end
end

this is useful to take advantages of broadcasting in the sense that you don't need to generate a dense meshgrid array but instead several "vectors" in orthogonal dimensions.

julia> g(x::Number, y::Number) = x+y
g (generic function with 1 method)

julia> g.(xvec(1:10, -5:1:5)...)
10×11 Array{Int64,2}:
 -4  -3  -2  -1  0   1   2   3   4   5   6
 -3  -2  -1   0  1   2   3   4   5   6   7
 -2  -1   0   1  2   3   4   5   6   7   8
 -1   0   1   2  3   4   5   6   7   8   9
  0   1   2   3  4   5   6   7   8   9  10
  1   2   3   4  5   6   7   8   9  10  11
  2   3   4   5  6   7   8   9  10  11  12
  3   4   5   6  7   8   9  10  11  12  13
  4   5   6   7  8   9  10  11  12  13  14
  5   6   7   8  9  10  11  12  13  14  15

julia> X = collect(0:0.001:1)

julia> @btime xvec($X);
  159.338 ns (8 allocations: 272 bytes)

julia> @btime xvec($X, $X);
  281.436 ns (15 allocations: 592 bytes)

julia> @btime xvec($X, $X, $X);
  390.059 ns (20 allocations: 848 bytes)

julia> @btime xvec($X, $X, $X, $X);
  510.880 ns (25 allocations: 1.17 KiB)

P.S. Iterators.product almost does the same thing except that I need to define a tuple version of g:

g(args::Tuple) = g(args...)
g.(Iterators.product(1:10, -5:1:5))

cc: @stevengj

[johnnychen94]

Sorry for the noise, actually the comprehension list is a faster version:

julia> g.(xvec(X, Y)...) == [g(x,y) for x in X, y in Y]
true

julia> @btime g.(xvec($X, $X)...);
  549.074 μs (19 allocations: 7.65 MiB)

julia> @btime [g(x,y) for x in X, y in Y];
  8.609 μs (5 allocations: 86.20 KiB)