# 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
I think Flatten[Transpose[{xx, yy, c}, {3, 1, 2}], 1] == Flatten[Map[Part, Transpose[{xx, yy, c}, {3, 1, 2}], {2}], 1]. –  expression Aug 19 '13 at 9:00
@ruebenko yes that will work, but it will not include the actual x and y coordinates. It will only display the z-cordinates and will use for the x,y coordinates the number of points. I wanted it to be like with Matlab's where the x,y axis has the actual coordinate ticks. –  Nasser Aug 19 '13 at 9:20
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

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... –  The Toad Aug 19 '13 at 4:34