## Hump-day Challenger - MATLAB Indexing

### Matt Fig (view profile)

on 30 Mar 2011
Latest activity Answered by Jose Sanchez

### Jose Sanchez (view profile)

on 23 Oct 2018
Accepted Answer by Sean de Wolski

### Sean de Wolski (view profile)

This challenger is very easy to state:
Given an array of random size and dimensions whose only nonzero element is the first element, find the locations of the nonzero elements of the array returned by a random call to the REPMAT function.
Now for the details. Provided below is a test code for use in testing your solution to the challenger. In the code an array A is created, with random size and number of dimensions held in the vector SA. Another vector RA is created to be used as the second argument to REPMAT. Then a call is made to find the locations of the non-zero elements of the repmatted array using:
I = find(repmat(A,RA));
The challenge is to obtain the same vector as I, following two rules:
1. Your function must not reproduce the repmatted array, or any array of equal or larger size, either implicitly or explicitly, in any form. What this means practically is that your function will succeed when REPMAT fails due to either an out of memory problem or a maximum variable size exceeded problem.
2. Your function should be faster than the call to FIND and REPMAT shown above.
For those who are indexing experts this may not be much of a challenge, but for those who are still learning this could be very challenging indeed. I think one could learn a lot by completing the challenge, even if it is difficult, so keep trying! Study the test function below to understand the problem better and see how your code will be used (input args, etc). Remember, when testing your solution run the test function more than once in a row to get good timing results. Good luck. And don't forget to vote up good answers! .
.
.
function [] = test_solutions()
% Use to test an Answer for the Hump-Day Challenger.
N = 300; % The number of loop iterations.
Tm = [0 0]; % To hold the timings.
bf = 0; % Flag to check if user code passed.
Z = 6; % If you have consistent memory problems, lower this...
for ii = 1:N
% First make our array, and decide how to repmat it.
SA = ceil(rand(1,ceil(rand*Z))*(Z+1)); % The size of A.
A = false(SA);
A(1) = true;
RA = ceil(rand(1,ceil(rand*Z))*(Z+1)); % The vector for REPMAT.
% Now compare speeds.
tic
I = find(repmat(A,RA));
Tm(1) = Tm(1) + toc;
tic
Imf = hdchallenger(SA,RA); % Your func should take only SA and RA.
Tm(2) = Tm(2) + toc;
if ~isequal(I(:),Imf(:))
disp(' Unequal results, please try again.')
bf = 1;
break
end
end
RT = Tm/min(Tm); % The relative timing.
if RT(2)==1 && ~bf
fprintf('You passed! You beat FIND by a factor of: %.1f!\n',RT(1))
elseif ~bf
fprintf('Your code returns good results, but is slow. Keep trying!\n')
end

Matt Fig

### Matt Fig (view profile)

on 25 May 2011
Thanks to Sean de and Andrew Newell both!

### Sean de Wolski (view profile)

on 30 Mar 2011

function idx = hdchallenger(SA,RA)
%%SCd 03/30/2011
%Updated
if isscalar(RA)
RA = [RA RA];
end
if isscalar(SA)
SA = [SA SA];
end
lenRA = length(RA);
lenSA = length(SA);
if lenRA<lenSA
RA(lenSA) = 1;
RA(RA==0) = 1;
end
%Engine
n = prod(RA);
didx = SA(1)*ones(n-1,1); %preallocate
for ii = 2:lenRA;
skip = n/prod(RA(ii:end));
if ii > lenSA; %Don't exceed SA!
break
end
mpart = zeros(1,ii);
mpart(ii) = -1; %remove present 1 dim from count
didx(skip:skip:end) = prod([RA(1:ii-1),(SA(1:ii)+mpart)])+didx(skip:skip:end);
end
idx = cumsum([1; didx]); %integrate!

Matt Fig

### Matt Fig (view profile)

on 30 Mar 2011
Nice Sean de! This is faster than Andrew's engine on my machine. You too are a master MATLAB indexer!
Andrew Newell

### Andrew Newell (view profile)

on 30 Mar 2011
It's faster on my machine too. I get killed by the repmat's in my code. ### Andrew Newell (view profile)

on 31 Mar 2011

I did a lot of tweaking to get rid of repmat, preallocate II, and replace sub2ind by more efficient code. Alas, it is still slower than @Sean de's offering.
function I = hdchallenger1(SA,RA)
% repmat(A,m) = repmat(A,m,m)
if length(RA)==1
RA = [RA RA];
end
% false(m) = false(m,m)
if length(SA)==1
SA = [SA SA];
end
if length(RA)>length(SA)
SA = [SA ones(1,length(RA)-length(SA))];
else
RA = [RA ones(1,length(SA)-length(RA))];
end
% Calculate subscripts for ones
II = zeros(prod(RA),length(RA));
m1 = [1 cumprod(RA)];
m2 = m1(end)./m1;
for j=1:length(RA)
JJ = 1 + cumsum([zeros(m1(j),1) SA(j)*ones(m1(j),RA(j)-1)],2);
JJ = JJ(:);
JJ = JJ(:,ones(1,m2(j+1)));
II(:,j) = JJ(:);
end
% Convert to linear index (cf. sub2ind)
siz = RA.*SA;
k = [1 cumprod(siz(1:end-1))];
I = 1+ (II-1)*k';

Sean de Wolski

### Sean de Wolski (view profile)

on 31 Mar 2011
Would there be a way to feed your logic into NDGRID to have it compute the various combinations? Getting rid of the NUM2CELL conversion probably helped quite a bit too!
Andrew Newell

### Andrew Newell (view profile)

on 31 Mar 2011
Yes, getting rid of NUM2CELL did help a lot (although the cumulative effect of all the changes was disappointing). Your NDGRID idea is intriguing, but I'd better get back to my real work! You are welcome to give it a try. ### Andrew Newell (view profile)

on 30 Mar 2011

function I = hdchallenger(SA,RA)
% repmat(A,m) = repmat(A,m,m)
if length(RA)==1
RA = [RA RA];
end
% false(m) = false(m,m)
if length(SA)==1
SA = [SA SA];
end
if length(RA)>length(SA)
SA = [SA ones(1,length(RA)-length(SA))];
else
RA = [RA ones(1,length(SA)-length(RA))];
end
II = 1;
for j=1:length(RA)
JJ = (1:SA(j):(SA(j)*(RA(j)-1)+1));
JJ = repmat(JJ,size(II,1),1);
if j==1
II = JJ(:);
else
II = [repmat(II,RA(j),1) JJ(:)];
end
end
Ic = num2cell(II,1);
I = sub2ind(RA.*SA,Ic{:});
if SA(1)*RA(1)==1
I = I(:);
end

Sean de Wolski

### Sean de Wolski (view profile)

on 30 Mar 2011
What causes find to return a row vector? I can't seem to 'find' any pattern!
Matt Fig

### Matt Fig (view profile)

on 30 Mar 2011
FIND returns a row vector when the search array is a row vector or a scalar. I will change the test to:
~isequal(I(:),Imf(:))
Andrew Newell

### Andrew Newell (view profile)

on 30 Mar 2011
Thank you! That was driving me nuts! ### Andrew Newell (view profile)

on 31 Mar 2011

I tried to finesse this problem using the function sptensor from the Tensor Toolbox, but alas - repmat doesn't work properly because sptensor doesn't use matrix indices the way that ordinary multidimensional arrays do (see my question on Indexing arrays with matrices).

### Jose Sanchez (view profile)

on 23 Oct 2018

I left below my solution based on bsxfun:
function ind = hdchallenger(SA,RA)
if length(SA) == 1
SA = [SA SA];
end
if length(RA) == 1
RA = [RA RA];
end
ind = 1;
SA = [SA ones(1,max(0,length(RA)-length(SA)+1))];
ne = SA(1);
for it = 1:length(RA)
ind = bsxfun(@plus, ind, (0:RA(it)-1)*ne);
ind = ind(:);
ne = ne*RA(it)*SA(it+1);
end