# vectorize a for loop

8 visualizzazioni (ultimi 30 giorni)
Sam il 26 Ott 2012
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
-Sam
##### 1 CommentoMostra -1 commenti meno recentiNascondi -1 commenti meno recenti
José-Luis il 26 Ott 2012
Loops are not always evil, and in this case probably even faster than the alternatives.

Accedi per commentare.

### Risposta accettata

Sean de Wolski il 26 Ott 2012
Modificato: Sean de Wolski il 26 Ott 2012
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 :) .
##### 6 CommentiMostra 4 commenti meno recentiNascondi 4 commenti meno recenti
Sean de Wolski il 26 Ott 2012
Modificato: Sean de Wolski il 26 Ott 2012
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 il 27 Ott 2012
@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.

Accedi per commentare.

### Più risposte (5)

Matt J il 26 Ott 2012
Modificato: Matt J il 26 Ott 2012
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);
##### 8 CommentiMostra 6 commenti meno recentiNascondi 6 commenti meno recenti
Sean de Wolski il 26 Ott 2012
@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 il 26 Ott 2012
Yes, my tes is'nt adequate, like you said, I must make a test in a loop.

Accedi per commentare.

Azzi Abdelmalek il 26 Ott 2012
eval(sprintf('isPeak(start_peak_hour(%d):end_peak_hour(%d)) = 1,',[1:n_peak;1:n_peak]))
##### 3 CommentiMostra 1 commento meno recenteNascondi 1 commento meno recente
Azzi Abdelmalek il 26 Ott 2012
If you want to avoid eval
isPeak(cell2mat(arrayfun(@(x,y) x:y, start_peak_hour,end_peak_hour,'un',0)))=1
Jan il 26 Ott 2012
And you should avoid EVAL...

Accedi per commentare.

Jan il 26 Ott 2012
Modificato: Jan il 27 Ott 2012
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);
[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);
##### 3 CommentiMostra 1 commento meno recenteNascondi 1 commento meno recente
Jan il 26 Ott 2012
Modificato: Jan il 26 Ott 2012
At least I've considered the last element...
Matt J il 26 Ott 2012
At least I've considered the last element...
True!

Accedi per commentare.

Chris A il 26 Ott 2012
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'),
##### 0 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

Jan il 27 Ott 2012
Modificato: Jan il 27 Ott 2012
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
##### 0 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

### Categorie

Scopri di più su Introduction to Installation and Licensing in Help Center e File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by