1

I have a grid that I want to store x and y coordinates for in 1D arrays. It is orthogonal so that x values are the same at each y in the grid, and vice versa. I've recently tried to vectorize the procedure and am using a repmat-reshape combination to achieve this. I heard that repmat is not the best, and wondered if there are faster alternatives to this same procedure.

For example, if the grid contains 3 nodes in the x direction (nX = 3) and 4 in the y direction (nY = 4), and their "template" values are as follows:

xTemplate =   [0    0.5000    1.0000]
yTemplate =   [0    0.3333    0.6667    1.0000]

Then the values for the entire grid would look like:

xArray =
     0    0.5000    1.0000
     0    0.5000    1.0000
     0    0.5000    1.0000
     0    0.5000    1.0000

and

yArray =
     0         0         0
     0.3333    0.3333    0.3333
     0.6667    0.6667    0.6667
     1.0000    1.0000    1.0000

Ultimately, I want them to be contained in 1D arrays, where the values are stored in order from 1 to nX*nY, in other words where

iVec = nX*(jArray - 1) + iArray

So that they look as follows:

[xVec yVec] = 
     0         0
0.5000         0
1.0000         0
     0    0.3333
0.5000    0.3333
1.0000    0.3333
     0    0.6667
0.5000    0.6667
1.0000    0.6667
     0    1.0000
0.5000    1.0000
1.0000    1.0000

Currently, I have the following code to achieve this, but am wondering if there are other ways to increase the speed here:

nX = 3;
nY = 4;

xTemplate = linspace(0,1,nX); 
xArray    = repmat(xTemplate, nY, 1);
xVec      = reshape(xArray', [nX*nY 1]);

yTemplate = linspace(0,1,nY); 
yArray    = repmat(yTemplate, nX, 1)';
yVec      = reshape(yArray' , [nX*nY 1]);

What are some alternatives to this that might increase speed? Is there a simpler approach that omits creating xArray and yArray which are only intermediate and not necessary? Thanks for any help.

teepee
  • 2,620
  • 2
  • 22
  • 47
  • 2
    The two-column matrix you ultimately want to obtain can be generated using any of the [answers here](https://stackoverflow.com/q/21895335/2586922) with `vectors = {yTemplate, xTemplate}` and applying `fliplr` to the result. (This is also valid for more than two vectors) – Luis Mendo Nov 08 '17 at 18:38
  • thanks for that. It seems as though `meshgrid` is faster though, would that be something I did wrong? – teepee Nov 08 '17 at 19:57
  • I would expect plain `meshgrid` to quite fast here. My linked answer essentially uses that (well, the equivalent `ndgrid`), plus some additional code. Are you timing with `tic`/ `toc` or with `timeit`? The latter is more accurate – Luis Mendo Nov 08 '17 at 20:03
  • I've been using tic-toc. How does the other work better for timing, ooc? – teepee Nov 08 '17 at 21:27
  • [Here](https://stackoverflow.com/a/18955502/2586922)'s a nice explanation – Luis Mendo Nov 08 '17 at 22:29

3 Answers3

3

You can easily build xVec using repmat by just reshaping xTemplate to a column vector first using the colon operator (or a nonconjugate transpose .'):

xVec = repmat(xTemplate(:), [nY 1]);

And yVec can be built just as easily using repelem (present since R2015a):

yVec = repelem(yTemplate(:), nX);

If you're using an older version of MATLAB, this post gives a number of alternatives to repelem.

Timing:

I was curious, so I decided to time all the solutions thus far (with minor changes to each to make them more equivalent tests). Here are some typical timing results (using timeit, MATLAB R2016b, Windows 7), with the test code below:

>> grid_timing(3, 4)
teepee:  7.09515e-06  % Longest
gnovice: 4.66416e-06  % Shortest
Aero:    4.90961e-06
Dev-iL:  6.23629e-06

>> grid_timing(100, 100)
teepee:  4.41746e-05  % Longest
gnovice: 1.65685e-05  % Shortest
Aero:    3.00774e-05
Dev-iL:  2.81186e-05

The differences were often very small, even for larger grids, but my solution tends to be slightly faster.

Test code:

function grid_timing(nX, nY)
  [xt, yt] = teepee_grid(nX, nY);
  [xg, yg] = gnovice_grid(nX, nY);
  [xA, yA] = Aero_grid(nX, nY);
  [xD, yD] = DeviL_grid(nX, nY);
  if ~isequal(xt, xg, xA, xD) || ...
     ~isequal(yt, yg, yA, yD)
    disp('Solutions don''t match!');
  end
  fprintf('teepee:  %g\n', timeit(@() teepee_grid(nX, nY), 2));
  fprintf('gnovice: %g\n', timeit(@() gnovice_grid(nX, nY), 2));
  fprintf('Aero:    %g\n', timeit(@() Aero_grid(nX, nY), 2));
  fprintf('Dev-iL:  %g\n', timeit(@() DeviL_grid(nX, nY), 2));
end
function [xVec, yVec] = gnovice_grid(nX, nY)
  xVec = repmat(linspace(0,1,nX).', [nY 1]);
  yVec = repelem(linspace(0,1,nY).', nX);
end
function [xVec, yVec] = Aero_grid(nX, nY)
  yVec = linspace(0,1,nY).*ones(nX,1);
  xVec = linspace(0,1,nX).'.*ones(1,nY);
  yVec = yVec(:);
  xVec = xVec(:);
end
function [xVec, yVec] = DeviL_grid(nX, nY)
  [YY, XX] = meshgrid(linspace(0,1,nY), linspace(0,1,nX));
  xVec = XX(:);
  yVec = YY(:);
end
function [xVec, yVec] = teepee_grid(nX, nY)
  xArray    = repmat(linspace(0,1,nX), nY, 1);
  xVec      = reshape(xArray', [nX*nY 1]);
  yArray    = repmat(linspace(0,1,nY), nX, 1)';
  yVec      = reshape(yArray' , [nX*nY 1]);
end
gnovice
  • 125,304
  • 15
  • 256
  • 359
  • Thanks for that help, that is a really concise way to write it. However it doesn't seem to be faster so I'm wondering if there's any other tricks to add perhaps? – teepee Nov 08 '17 at 17:35
  • @teepee: How much time is this actually taking? I can't imagine this takes long enough to be a bottleneck in your code. Trying to speed this up seems like a micro-optimization. – gnovice Nov 08 '17 at 17:39
  • It has a larger effect when creating grids of thousands of nodes essentially, but I think you're right that it shouldn't be a bottleneck. – teepee Nov 08 '17 at 17:45
  • That is sweet. Very nice look at it there. Thanks for doing that. Also, I noticed that using `reshape` in my function caused it to crash compared with your (:) operations for very large grids. How does the memory treatment work differently in those cases, ooc? – teepee Nov 08 '17 at 21:30
  • @teepee: The memory issue may simply be due to the fact that you are storing some intermediate values (i.e. `xArray` and `yArray`) which may be taking up more memory, and it just so happens that MATLAB runs out of memory during the reshape operation. – gnovice Nov 08 '17 at 21:40
3

To make the code shorter, you can use meshgrid for this:

nX = 3;
nY = 4;

[YY, XX] = meshgrid(linspace(0,1,nX), linspace(0,1,nY));
xVec = XX(:);
yVec = YY(:);

Though this probably doesn't improve performance...

Dev-iL
  • 23,742
  • 7
  • 57
  • 99
  • This is actually the fastest one yet according to some tests – teepee Nov 08 '17 at 18:49
  • You need to swap your input arguments to `meshgrid`, too. – gnovice Nov 08 '17 at 19:10
  • I noticed that as well. It's a bit confusing for readability to have `[Y,X] = meshgrid(y,x)`, but it's concise and fast enough that I might stick with it/ – teepee Nov 08 '17 at 19:57
  • @teepee: How did you do the timing? Using [`timeit`](https://www.mathworks.com/help/matlab/ref/timeit.html) is generally more accurate than [`tic`](https://www.mathworks.com/help/matlab/ref/tic.html)/[`toc`](https://www.mathworks.com/help/matlab/ref/toc.html). – gnovice Nov 08 '17 at 20:30
  • I did use tic-toc. Oops. How far off is tic-toc usually? – teepee Nov 08 '17 at 21:28
  • @teepee: Some comparisons [here](https://www.mathworks.com/help/matlab/matlab_prog/measure-performance-of-your-program.html). `timeit` takes the median over multiple runs and accounts for first-time costs. `tic/toc` also isn't recommended for fast code (less than 0.1 sec). – gnovice Nov 08 '17 at 21:54
2

Here is an alternative using implicit expansion in newer version of Matlab:

yVec = yTemplate .* ones(nX,1);
xVec = xTemplate' .* ones(1,nY);

yVec = yVec(:);
xVec = xVec(:);
Aero Engy
  • 3,588
  • 1
  • 16
  • 27