Meshgrid function

[mauri...@gmail.com]

Is the meshgrid function (similar to MATLAB) available in Julia?

I'm trying to plot a surface using 3 (1D) vectors and I'm not sure yet how can I deal with this problem
without using this function.

Thank you!

[Isaiah Norton]

See:

https://groups.google.com/forum/#!searchin/julia-users/meshgrid|sort:date/julia-users/FuKK7zjncN8/dpWjRImMBgAJ

[mauri...@gmail.com]

Thank you for your note,

But the problem is that broadcasting will not work in my case, once I already have the three vectors that I want to plot instead of
two entries and one image (as in the example). I believe that I need something to interpolate my vectors in order to generate the surface.

[Isaiah Norton]

I think most pure-Julia interpolation work has converged to the Interpolations package: https://github.com/tlycken/Interpolations.jl

(several alternatives are also listed on that page)

[Cedric St-Jean]

I used the meshgrid in Julia's examples folder: https://github.com/JuliaLang/julia/blob/master/examples/ndgrid.jl

[Tim Holy]

Try
Ainterp = [itp[x,y] for x in X, y in Y]

It's best to do this inside a function, because then the element type should
be inferred, and it will be far more efficient.

Best,
--Tim

On Sunday, March 06, 2016 05:22:17 PM mauri...@gmail.com wrote:
> Thank you for your note,
>
> But the problem is that broadcasting will not work in my case, once I
> already have the three vectors that I want to plot instead of
> two entries and one image (as in the example). I believe that I need
> something to interpolate my vectors in order to generate the surface.
>
> On Sunday, March 6, 2016 at 7:17:38 PM UTC-6, Isaiah wrote:
> > See:
> >
> > https://groups.google.com/forum/#!searchin/julia-users/meshgrid|sort:date/
> > julia-users/FuKK7zjncN8/dpWjRImMBgAJ>

[Christoph Ortner]

Have a look at 

  http://matlabcompat.github.io

If it doesn't have a mesh grid function, then maybe mesh grid from the examples folder could be added to this package. 

Christoph

[Tomas Lycken]

The thing with meshgrid, is that it’s terribly inefficient compared to the equivalent idioms in Julia. Instead of implementing it and keeping with inefficient habits, embrace the fact that there are better ways to do things :)

Since this question is recurring, I think the problem might rather be documentation. As a user looking for a meshgrid function, where did you look? What did you look for, more specifically? Would it have helped to have some examples along the lines of “when you want to do this, using meshgrid, in MATLAB, here’s how to do the same thing without meshgrid in Julia”? I could turn the following writeup into a FAQ entry or something, but it would be nice to know where to place it to actually have it make a difference.

Making do without meshgrid (and getting out ahead!)

For starters, if what you are after is a straight-up port, meshgrid is very easy to implement for yourself, so that you don’t have to think about the rest. Mind you, this is not the most performant way to do this!

meshgrid(xs, ys) = [xs[i] for i in 1:length(xs), j in length(ys)], [ys[j] for i in 1:length(xs), j in 1:length(ys)]

Try it out with X, Y = meshgrid(1:15, 1:8). (If you’re not happy with the orientation of the axes, just switch i and j in the meshgrid implementation.)

However, most of the uses of meshgrid in MATLAB is for stuff where you want to do a bunch of calculations for each point on a grid, and where each calculation is dependent on the coordinates of just one point. Usually something along the lines of

X, Y = meshgrid(xs, ys)
Z = sin(X) + 2 * Y.^2 ./ cos(X + Y)

In Julia, this is better done with a list comprehension altogether, eliminating the need for allocating X and Y:

Z = [sin(x) + 2 * y^2 / cos(x+y) for x in xs, y in ys]

If you think it’s now become too difficult to read what the actual transformation is, just refactor it out into a function:

foo(x,y) = sin(x) + 2 * y^2 / cos(x+y)
Z = [foo(x,y) for x in xs, y in ys]

If you do this a lot, you might want to avoid writing the list comprehension altogether; that’s easy too. Just create a function that does that for you. With a little clever interface design, you’ll see that this is even more useful than what you could do with meshgrid.

inflate(f, xs, ys) = [f(x,y) for x in xs, y in ys]

# call version 1: with a named function
inflate(foo, xs, ys) 

# call version 2: anonymous function syntax
inflate((x,y) -> sin(x) + 2 * y^2 / cos(x,y), xs, ys)

# call version 3: with a do block
inflate(xs, ys) do x, y
    sin(x) + 2 * y^2 / cos(x,y)
end

Note that the last version of this call is very flexible; in fact, it’s just as flexible as having a named function. For example, you can do everything in multiple steps, allowing for intermediate results that don’t have to allocate large arrays (as they would have, had you used meshgrid), making it easier to be explicit about whatever algorithm you’re using.

However, the two last versions of these calls are not the most performant as of yet, since they’re using anonymous functions. This is changing, though, and in 0.5 this won’t be an issue anymore. Until then, for large arrays I’d use the first call version, with a named function, to squeeze as much as possible out of this.

I did a small benchmark of this using a straight list comprehension (i.e. without the inflate function) vs using meshgrid. Here are the results:

julia> foo(x,y) = sin(x) + 2y.^2 ./ cos(x+y);

julia> function bench_meshgrid(N)
           X, Y = meshgrid(1:N, 1:N)
           foo(X, Y)
       end;

julia> function bench_comprehension(N)
           [foo(x,y) for x in 1:N, y in 1:N]
       end;

# warmup omitted

julia> @time bench_comprehension(1_000);
  0.033284 seconds (6 allocations: 7.630 MB)

julia> @time bench_meshgrid(1_000);
  0.071993 seconds (70 allocations: 68.666 MB, 5.12% gc time)

julia> @time bench_comprehension(10_000);
  3.547387 seconds (6 allocations: 762.940 MB, 0.06% gc time)

julia> @time bench_meshgrid(10_000);
ERROR: OutOfMemoryError()
 in .* at arraymath.jl:118
 in foo at none:1
 in bench_meshgrid at none:3

Note the memory usage: for the list comprehension, it’s entirely dominated by the result (10_000^2 * 8 / 1024^2 ~= 762.939 MB), while meshgrid eats almost ten times as much memory (the exact ratio will vary very much depending on the function foo and how many intermediate arrays it creates). For large inputs, meshgrid doesn't even make it.

// T

[Eric Forgy]

That's awesome Tomas. Thanks for the lesson :D

[ele...@gmail.com]

Neat.

Wherever it goes it needs a pointer from http://docs.julialang.org/en/latest/manual/noteworthy-differences/#noteworthy-differences-from-matlab but its a bit long to actually go there.

[mauri...@gmail.com]

Thank you very much, it was really helpful!

[Christoph Ortner]

Nice list of alternatives. 

But it doesn't change the fact that meshgrid is very easy to use, and very easy to read. New idioms should not be forced on developers. I know in what ways mesh grid is inefficient and when that matters then I won't use it. But more often than not I just want a quick have for a simple problem, and then I will use meshgrid over list comprehensions, etc. any time.

Christoph

[Tobias Knopp]

I also don't see the major issue with the inefficiency of meshgrid.

We do provide functions for evaluating

c = a * b

with a, b being vectors and we know that it is not efficient to do so (due to the temporary vector introduced)

Tobi

[Tomas Lycken]

Temporary arrays for other things, like C = A * B (or, more relevant, C = A .* B), is an optimization issue actively being worked on. But it’s a hard one to solve, unlike avoiding meshgrid. Furthermore; if you avoid using meshgrid, you’ll avoid falling on this performance bottleneck too, as you will probably be able to avoid this type of array operation entirely.

Since building a meshgrid function yourself is so easy (and will be just as fast and as good as if it was included in Base), but using it is potentially harmful to performance, I think there are good reasons to keep it out of Base. If you really want to use it, it’s a one-liner in .juliarc away, and as long as it’s not in Base there’s a good chance that beginner Julians will pick up more performant idioms.

I'll look into posting this somewhere less transient, to make it available for linking from the Julia vs MATLAB FAQ.

// T

[Christoph Ortner]

" and as long as it’s not in Base there’s a good chance that beginner Julians will pick up more performant idioms."

Hi Tomas,

I think this is a big mistake so many Julia developers (not all) are making: Julia should not just be great for performance of code but also for performance of the programmer! Meshgrid just makes for very simple and readable code. 

Christoph

[J Luis]

I think this is a big mistake so many Julia developers (not all) are making:

I agree with this. For all those in future that are not necessarily programmers (and Julia hopes to attract many of those, I believe) and will need a meshgrid, the time will take them to find out the alternative can hardly be seen a efficiency factor from the language.
That said, meshgrid is one of my favorite hates in Matlab  so big is the resources waste it carries with, but when I need it I like it to be there.

[STAR0SS]

Julia devs are working hard trying to make a language than is consistent and elegant, they don't want to add functions in Base only because they are convenient for people coming from a particular background when there's a more "Julian" way of doing it. Plus if you start adding matlab's function, why stop there? Add also R and python ones, you see it's not really a practical solution.

The solution is to put all these in the matlabcompat package, so that people that are familiar with matlab can have all the familiar tools they need (meshgrid, interp1, polyfit, ...) without having to find, install and read the docs of 10 different packages. Since these functions are relatively simple, I'd encourage you to contribute to that package.

[Christoph Ortner]

On Wednesday, 9 March 2016 19:47:11 UTC, STAR0SS wrote:

Julia devs are working hard trying to make a language than is consistent and elegant, they don't want to add functions in Base only because they are convenient for people coming from a particular background when there's a more "Julian" way of doing it. Plus if you start adding matlab's function, why stop there? Add also R and python ones, you see it's not really a practical solution.

The solution is to put all these in the matlabcompat package, so that people that are familiar with matlab can have all the familiar tools they need (meshgrid, interp1, polyfit, ...) without having to find, install and read the docs of 10 different packages. Since these functions are relatively simple, I'd encourage you to contribute to that package.

This is precisely what I suggested a few posts ago.

Christoph