Matlab—replicate columns of a matrix unequal numbers of times

expando.m

function [] = expando()
  widthP = 1000;
  maxRep = 40;
  reps = randi(maxRep, 1, widthP) - 1;

  P = randn(10, widthP);
  widthM = sum(reps);

  M = method1(reps, P);

  M2 = method2(reps, P);
  assert(all(M2(:) == M(:)))

  [M3, i3] = method3(reps, P);
  assert(all(M3(:) == M(:)))

  disp('All asserts pass!')

  t1 = timeit(@() method1(reps, P));
  t2 = timeit(@() method2(reps, P));
  t3 = timeit(@() method3(reps, P));
  relativeTimes = [t1 t2 t3] / t1
end

function M = method1(reps, P)
  widthP = length(reps);
  M = [];
  for colIdx = 1 : length(reps)
    M = [M repmat(P(:, colIdx), 1, reps(colIdx))];
  end
end

function M = method2(reps, P)
  widthP = length(reps);
  c = arrayfun(@(n, i) i*ones(n,1), reps, 1:widthP, 'un', 0);

  i = cell2mat(c(:));
  M = P(:, i);
end

function [M3, i] = method3(reps, P)
  nzidx = reps~=0;
  nzreps = reps(nzidx);

  % reps != 0 algorithm as above, except using nzreps instead of reps
  widthM = sum(reps);
  z = zeros(widthM, 1);
  j = [1 1+cumsum(nzreps(1:end-1))]; % pixie dust, not much
  z(j) = 1;
  i = cumsum(z);

  % adjust i
  n = 1:numel(reps);
  f = n(nzidx);
  i2 = f(i);
  % here's another way to do that with `find` (which can be slow, and indeed is)
  if 0
    f = find(nzidx);
    i3 = f(i);
    assert(all(i2==i3))
  end

  % All done. Grab the appropriate columns of the input
  M3 = P(:, i2);
end
% relativeTimes =
%            1       1.1411     0.076827
%            1       1.2945     0.087953
%            1       1.1663      0.10526