# Simulate MATLAB's meshgrid function

Here is MATLAB's meshgrid. I came up with these implementations in Mathematica:

{x, y} = With[{n = 10}, Table[#, {i, n}, {j, n}]] & /@ {j, i};


and

{x, y, z} = Array[List, {5, 5, 5}][[All, All, #]] & /@ {3, 2, 1}


but it seems that it's not fast enough. The following example is about three times slower than MATLAB.

(* Mathematica *)
(
{a, b} = With[{n = 3000}, Table[#, {i, n}, {j, n}]] & /@ {j, i};
c = a^2 + b^2 // N // Sqrt;
Compile[{}, Position[Unitize@FractionalPart@c, 0]][] // Length
) // AbsoluteTiming

%MATLAB
tic;
[a,b]=meshgrid(1:3000);
c=sqrt(a.^2+b.^2);
idx=find(c==round(c));
length(idx)
toc;


Can you recommend an efficient method?

Updated

(
{a, b} = {#, Transpose@#} &@ConstantArray[Range@3000, {3000}];
c = Sqrt@N[a*a + b*b];
SparseArray[Unitize@FractionalPart@c~BitXor~1]["NonzeroPositions"] // Length
) // AbsoluteTiming

• For those of us who do not use Matlab can you describe what the meshgrid function does ...rather than forcing us to infer it by examining your code. – Mike Honeychurch Aug 19 '13 at 4:41

(Update, added more points, and more timings)

Using MATLAB's help standard example for meshgrid:

### Mathematica implementation

meshgrid[x_List, y_List]:={ConstantArray[x,Length[x]],Transpose@ConstantArray[y,Length[y]]}

{xx, yy} = meshgrid[Range[-2, 2, .1], Range[-4, 4, .2]];
c = xx*Exp[-xx^2 - yy^2];

pts = Flatten[{xx, yy, c}, {2, 3}];
ListPlot3D[pts, PlotRange -> All, AxesLabel -> Automatic,
ImagePadding -> 20, Mesh -> 35, InterpolationOrder -> 2,
ColorFunction -> "Rainbow", Boxed -> False]


### MATLAB

    [X,Y] = meshgrid(-2:.1:2, -4:.2:4);
Z = X .* exp(-X.^2 - Y.^2);
surf(X,Y,Z)


MATLAB timing was done using this template code by changing the grid spacing, as an example for the code that generate the above plot:

tic;  [X,Y] = meshgrid(-2:.2:2, -4:.4:4); toc


Mathematica was done on this template

Timing[
{xx, yy} = meshgrid[Range[-2, 2, .2], Range[-4, 4, .4]];
][[1]]


## Table of timings

The timing were done for different grid sizes as in (.2,.4) pairs in the above example and using Mathematica 9.01 and MATLAB 2012a, all on Windows 7. I could not do smaller grids than shown in the table below, else I ran out of memory. All times are in seconds. Same PC. 16 GB RAM, Intel Core i7

Mathematica's above implementation of meshgrid using ConstantArray seems to be a little slower as the number of grid points increased. A faster implementation should be investigated.

• I think ListPlot3D[Transpose[c],...] is enough. – user21 Aug 19 '13 at 8:52
• You could also use pts = Flatten[{xx, yy, c}, {2, 3}] – Simon Woods Aug 20 '13 at 10:34
• @SimonWoods thanks! that is even simpler. Just updated. I must say I find Matlab little easier when it comes to this part. Matlab wants each coordinates in one separate vector, so one passes 3 separate vectors, while Mathematica wants the coordinates in a matrix (each row has each {x,y,z} coordinate). I never understood why can't there be 2 API's for this in Mathematica. Where the user can choose to pass the coordinates as {X...},{Y...},{Z..} vs {{X,Y,Z}...}, and let the user choose the simple one depending. What is so hard about having two different interfaces like this for the command? – Nasser Aug 20 '13 at 14:56
• You can also use ListPlot3D[Transpose@c, DataRange -> {{-4, 4}, {-2, 2}}, PlotRange -> All] – chyanog Sep 2 '13 at 2:29
• At large n I think it's more efficient to build the list via Map than to Transpose it. i.e. use Map[ConstantArray[#, Length[y]] &, y] instead of Transpose@ConstantArray[y, Length@y]. It's not much better, but it does give a second or so of improvement when you have 15000 mesh points in both x and y – b3m2a1 Feb 3 '18 at 0:11

Concise

Another way to define meshgrid() in Mathematica is:

meshgrid[xgrid_List, ygrid_List] := Transpose[Outer[List, xgrid, ygrid], {3, 2, 1}]


It's definitely not winning the speed competition (about 10x slower than rm -rf's solution), but speaks more for Mathematica's language, which allows neat and concise definitions like this.

Speed

Regarding speed, in your particular example you can do the squaring before the repetition via ConstantArray, loose the superfluous N in the computation of c and use Count instead of Position[...]\\Length to save computation time.

({a, b} = {#, Transpose@#} &@ConstantArray[N@Range@3000, {3000}];
c = a^2 + b^2 // N // Sqrt;
Compile[{}, Position[Unitize@FractionalPart@c, 0]][] //Length) // AbsoluteTiming
(* {0.698965,7388} *)

({asq, bsq} = {#, Transpose@#} &@ ConstantArray[N[Range[3000]]^2, {3000}];
c = Sqrt[asq + bsq];
Compile[{}, Count[Unitize@FractionalPart@c, 0, 2]][]) // AbsoluteTiming
(* {0.396491,7388} *)

• (+1) I contemplated posting the same thing using Outer (which IMO is the clearest way to think about tensored grids) but I noticed that this turned into a speed chase and got discouraged. – gpap Aug 20 '13 at 16:42

meshgrid()'s functionality is easily handled by ConstantArray. For the example you provided (i.e. meshgrid called with a single vector),

{a, b} = {#, Transpose@#}&@ConstantArray[N@Range@3000, {3000}]; // AbsoluteTiming
(* {0.068030, Null} *)


You can easily extend it to handle the two argument case by calling ConstantArray for each list (x and y) and using the length of the other list as the number of repetitions, which I'll leave for you to implement.

• This brings down the time from ~0.9s to ~0.5s on my system. Still slower than MATLAB... I think an explicit Position finder (i.e., using Do loops) that is compiled to C in Mathematica might speed things up a bit more... – rm -rf Aug 19 '13 at 4:34

The accepted answer excellent and comprehensive, but it has a very fine/little mistake\bug. I've tried to fix it, but it has to be accepted due to my current reputation(=101) on this site(well, I'm not very active here), and mistakenly it was not. The mistake is in the correct treatment of lists of different sizes. I'll explain in more details...

consider two implementations, the first (meshgrid1) is copied from current version (as of 2/2/2018) of the accepted answer, the second (meshgrid2) is a correctly reproduces matlab's meshgrid. The little difference is in the argument of Length.

meshgrid1[x_List, y_List]:={ConstantArray[x,Length[x]],Transpose@ConstantArray[y,Length[y]]}

meshgrid2[x_List, y_List]:={ConstantArray[x,Length[y]],Transpose@ConstantArray[y,Length[x]]}


now use it in a slightly different example then in the accepted answer (the ylist here is not the same size as xlist)

xlist = Range[-2, 2, .1];
ylist = Range[-4, 4, .1];
Length[xlist]
Length[ylist]

Out: 41
Out: 81


and see the dimensions of the result

{xx1, yy1} = meshgrid1[xlist, ylist];
Dimensions[xx1]
Dimensions[yy1]

Out: {41, 41}
Out: {81, 81}


and

{xx2, yy2} = meshgrid2[xlist, ylist];
Dimensions[xx2]
Dimensions[yy2]

Out: {81, 41}
Out: {81, 41}


now the same thing in matlab - decide by yourself which one is correct or make more sense

>> xlist = -2:.1:2;
>> ylist = -4:.1:4;
>> [xx,yy] = meshgrid(xlist, ylist);

>> size(xlist)

ans =

1    41

>> size(ylist)

ans =

1    81
size(xx)

ans =

81    41

>> size(yy)

ans =

81    41

• Good job on debugging. The accepted answer probably got confused since the two starting vectors happened to be the same length, so the wrong lengths being supplied was not noticed. Had a rectangular plotting region been used like in your case, the error would have been spotted more quickly. – J. M. will be back soon Mar 21 '18 at 4:59