Simulate MATLAB's meshgrid function

Question

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

Answer 1

(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.

Answer 2

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} *)

Answer 3

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.

Answer 4

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