vectorize a for loop

I am finding a optimal way of vectorizing the following piece of code (of course the real data is much greater than this):

start_peak_hour = [3 10 20 50];
end_peak_hour   = [5 14 27 70]; % always same size as start_peak_hour
n_peak = numel(start_peak_hour);

isPeak = zeros(1,80); 

for j=1:n_peak,
    isPeak(start_peak_hour(j):end_peak_hour(j)) = 1;
end

Your two cents are appreciated

-Sam

[José-Luis]

Loops are not always evil, and in this case probably even faster than the alternatives.

Answer 1

[Sean de Wolski]

I highly doubt you'll find anything faster than that well-written for-loop. The JIT will pick that up and have wonders with it.

eval is evil and will be orders of magnitude slower and hideous/hard to follow.

start_peak_hour = [3 10 20 50];
end_peak_hour   = [5 14 27 70]; % always same size as start_peak_hour
n_peak = numel(start_peak_hour);

t1 = 0;
t2 = 0;
for tt = 1:50 %50x

tic
isPeak = zeros(1,80); 
for j=1:n_peak,
  isPeak(start_peak_hour(j):end_peak_hour(j)) = 1;
end
t1 = t1+toc;
tic
eval(sprintf('isPeak(start_peak_hour(%d):end_peak_hour(%d)) = 1,',[1:n_peak;1:n_peak]))
t2=t2+toc;
end
[t1 t2] % 0.0010    0.0681
t1/t2

On my system the FOR-loop is 70x faster (for this small problem), it'll get better with bigger ones. And it doesn't use eval or any other impossible to understand syntax.

My $0.05 - a full nickel :) .

[Sam]

thanks Sean, this helps alot. I am supprised that EVAL is much slower than for loop

[Sean de Wolski]

eval can not be accelerated by the JIT, thus it needs to be re-figured out each time and do all sorts of kludgy conversions. In almost five years of using MATLAB I've never found one good use for it.

[Sean de Wolski]

My answer here may interest you as well:

http://www.mathworks.com/matlabcentral/answers/40977#answer_50587

[Azzi Abdelmalek]

Sometimes we don't need speed, like in this example.

[Sean de Wolski]

Well I guess it depends on how we define Sam's "optimal way"

Some possible definitions:

• easy to read: goes to for-loop
• fast: goes to for-loop
• easy to debug: goes to for-loop
• memory efficient: tie
• Uses JIT features: goes to for-loop

So unless you want to confuse people and write cruddy code - avoid eval

[Jan]

@Azzi: I prefer to write fast code for all cases. A useful piece of code is used in larger programs also and then the total speed is limited by the most time consuming parts. A not useful piece of code is deleted soon.

A severe problem is the dependency to the problem size: while implementing a quicksort seems to be efficient, it is a bad idea for lists of less than 10 elements. Calculating a SUM using multiple threads and all cores is reasonable for large arrays, it is a waste of time for tiny arrays.

So I respond: We always need "speed", but the "speed" of an algorithm is not well defined.

Answer 2

[Matt J]

Here's a way that avoid both EVAL and for-loops

isPeak = sparse(80,1,sum(end_peak_hour-start_peak_hour+1));

ispeak(start_peak_hour)=1;  
ispeak(end_peak_hour+1)=-1; 
ispeak=cumsum(ispeak);

[Azzi Abdelmalek]

I am not getting the result

[Matt J]

Works fine for me.

[Matt J]

Beter yet,

e=ones(1,n_peak);
nzmax=sum(end_peak_hour-start_peak_hour+1);
ispeak = sparse([start_peak_hour,end_peak_hour+1],1, [e, -e],80,1,nzmax);

ispeak=cumsum(ispeak);

[Sean de Wolski]

Adding this into the time test, the for-loop is still 5x faster!

start_peak_hour = [3 10 20 50];
end_peak_hour   = [5 14 27 70]; % always same size as start_peak_hour
n_peak = numel(start_peak_hour);

t1 = 0;
t2 = 0;
t3 = 0;
for tt = 1:50 %50x

tic
isPeak = zeros(1,80); 
for j=1:n_peak,
    isPeak(start_peak_hour(j):end_peak_hour(j)) = 1;
end
t1 = t1+toc;
tic
eval(sprintf('isPeak(start_peak_hour(%d):end_peak_hour(%d)) = 1,',[1:n_peak;1:n_peak]))
t2=t2+toc;
tic
e=ones(1,n_peak);
nzmax=sum(end_peak_hour-start_peak_hour+1);
ispeak = sparse([start_peak_hour,end_peak_hour+1],1, [e, -e],80,1,nzmax);
ispeak=cumsum(ispeak);
t3=t3+toc;
end
[t1 t2 t3] % 0.0005    0.0260    0.0023

[Azzi Abdelmalek]

tic
isPeak = zeros(1,80); 
isPeak(cell2mat(arrayfun(@(x,y)  x:y,  start_peak_hour,end_peak_hour,'un',0)))=1;
toc
tic
isPeak = zeros(1,80); 
for j=1:n_peak,
isPeak(start_peak_hour(j):end_peak_hour(j)) = 1;
end
toc

Elapsed time is 0.027788 seconds.
Elapsed time is 0.028316 seconds.

[Matt J]

Depends a bit on the data, and the density of the intervals. Here's a case where the cumsum approach outperforms the loop.

 N=1e7;
 start_peak_hour = 1:20:N;
 end_peak_hour   = start_peak_hour+5; % always same size as start_peak_hour

 n_peak = numel(start_peak_hour);
 ispeak=zeros(1,N);

 tic;  

    ispeak(start_peak_hour)=1;
    ispeak(end_peak_hour+1)=-1;
    ispeak=cumsum(ispeak);

 toc   %Elapsed time is 0.031681 seconds.

 tic

   for j=1:n_peak,
       isPeak(start_peak_hour(j):end_peak_hour(j)) = 1;
   end

 toc   %Elapsed time is 0.374747 seconds.

[Sean de Wolski]

@Matt, fair enough! That gets your method a 6x speedup over mine.

@Azzi, that isn't really a fair test since you're only calling it once so sinmple java things can throw the results way off and not get all of the JIT benefits:

N = 80;
start_peak_hour = [3 10 20 50];
end_peak_hour   = [5 14 27 70]; % always same size as start_peak_hour
n_peak = numel(start_peak_hour);

%  N=1e7;
%  start_peak_hour = 1:20:N;
%  end_peak_hour   = start_peak_hour+5; % always same size as start_peak_hourn_peak = numel(start_peak_hour);
%  n_peak = numel(start_peak_hour);
%  
t1 = 0;
t2 = 0;
t3 = 0;
for tt = 1:10 %50x

tic
isPeak=zeros(1,N);
for j=1:n_peak,
    isPeak(start_peak_hour(j):end_peak_hour(j)) = 1;
end
t1 = t1+toc;
tic
isPeak = zeros(1,80); 
isPeak(cell2mat(arrayfun(@(x,y)  x:y,  start_peak_hour,end_peak_hour,'un',0)))=1;
t3=t3+toc;
end
[t1 t3]

I See a 27x speedup with the for-loop. Remember, converting to and from cells is slow and array/cellfun is typically slower than a for-loop too.

[Azzi Abdelmalek]

Yes, my tes is'nt adequate, like you said, I must make a test in a loop.

Answer 3

[Azzi Abdelmalek]

 eval(sprintf('isPeak(start_peak_hour(%d):end_peak_hour(%d)) = 1,',[1:n_peak;1:n_peak]))

[Sam]

thanks Azzi

[Azzi Abdelmalek]

If you want to avoid eval

isPeak(cell2mat(arrayfun(@(x,y)  x:y,  start_peak_hour,end_peak_hour,'un',0)))=1

[Jan]

And you should avoid EVAL...

Answer 4

[Jan]

[EDITED] See Matt J's answer, which I had overseen during typing.

start_peak_hour = [3 10 20 50];
end_peak_hour   = [5 14 27 70]; % always same size as start_peak_hour

isPeak = zeros(1,80);
isPeak(start_peak_hour)   = 1;
if end_peak_hour(end) == 80
  end_peak_hour(end) = [];
end
isPeak(end_peak_hour + 1) = -1;

isPeak = cumsum(isPeak);

Please measure the timings with your original data set.

[EDITED] A new approach:

[EDITED 2]: Bugs fixed, now it compiles:

#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[]) {
mwSize i, n, len;
double *ini, *fin, *r, *q, *qf;

ini = mxGetPr(prhs[0]);
fin = mxGetPr(prhs[1]);
n   = mxGetNumberOfElements(prhs[0]);
len = (mwSize) mxGetScalar(prhs[2]);

plhs[0] = mxCreateDoubleMatrix(1, len, mxREAL);
r       = mxGetPr(plhs[0]) - 1;

for (i = 0; i < n; i++) {
  q  = r + (mwSize) ini[i];
  qf = r + (mwSize) fin[i];
  while(q <= qf) {
    *q++ = 1.0;
  }
}

return;
}

In a real-world program checks of the inputs are OBLIGATORY.

Compile it, call it as:

myPeakFiller(start_peak_hour,end_peak_hour,80);

[Matt J]

Ummm. That's the same as mine, I'm afraid. :-)

In any case, it's pretty clear that it can only beat a for-loop if there is a sufficient density of intervals. Sean and I also showed this by example in the stream of comments to my answer.

[Jan]

Of course, sorry, Matt J. I've open the thread some hours ago, had something to work (problems with the Windows Update on the laptop...), and posted the answer without re-loading the thread. Therefore I haven't seen your answer before and because I obviously prefer this method, I'm going to vote your answer now.

At least I've considered the last element...

[Matt J]

At least I've considered the last element...

True!

Answer 5

[Chris A]

start_peak_hour = [3 10 20 50]';
end_peak_hour   = [5 14 27 70]'; % always same size as start_peak_hour
peak_size = end_peak_hour - start_peak_hour + 1;
ncols=max(peak_size);
nrows=numel(peak_size);
tmpmat=kron(start_peak_hour,ones(1,ncols))+(cumsum(ones(nrows, ncols), 2)-1),
lbmat=kron(start_peak_hour-1, ones(1,ncols));
ubmat=kron(end_peak_hour+1, ones(1,ncols));
u=tmpmat((tmpmat>lbmat)&(tmpmat<ubmat));
isPeak = zeros(1,80); 
isPeak(u) = 1;
horzcat((1:80)',isPeak'),

Answer 6

[Jan]

List of solutions:

function isPeak = Loop_(a, b, len)
isPeak = zeros(1, len); 
for j = 1:length(a)
  isPeak(start_peak_hour(j):end_peak_hour(j)) = 1;
end

%
function isPeak = CumSum_(a, b, len)
isPeak = zeros(1,80);
isPeak(a) = 1;
if b(end) == 80
  b(end) = [];
end
isPeak(b + 1) = -1;
isPeak = cumsum(isPeak);

%
function isPeak = Array_(a, b, len)
isPeak(cell2mat(arrayfun(@(x,y)  x:y, a, b, 'un', 0))) = 1;

%
function isPeak = Eval_(a, b, len)
eval(sprintf('isPeak(a(%d):b(%d)) = 1;' , [1:n_peak;1:n_peak]));

%
function isPeak = CumSum2_(a, b, len)
e      = ones(1, n_peak);
nzmax  = sum(b - a + 1);
isPeak = cumsum(sparse([a, b+1],1, [e, -e], len, 1, nzmax));

%
% And the C-Mex I've posted before

A larger test data set:

start_peak_hour = [3 10 20 50];
end_peak_hour   = [5 14 27 70];
a = bsxfun(@plus, start_peak_hour', 0:100:100000);
a = a(:)';
b = bsxfun(@plus, end_peak_hour', 0:100:100000);
b = b(:)';
len = b(end) + 10;

Timings:

tic; for k = 1:1000, isPeak = FUNC(a, b, len); end; toc

Loop_:      4.75 sec
CumSum_:    1.19 sec
Array_:    37.08 sec
Eval_:   1224.83 sec
CumSum2_:   4.09 sec
Mex:        0.42 sec