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
1 个评论
José-Luis
2012-10-26
Loops are not always evil, and in this case probably even faster than the alternatives.
采纳的回答
Sean de Wolski
2012-10-26
编辑：Sean de Wolski
2012-10-26
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 个评论
Sam
2012-10-26
Sean de Wolski
2012-10-26
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
2012-10-26
My answer here may interest you as well:
Azzi Abdelmalek
2012-10-26
Sometimes we don't need speed, like in this example.
Sean de Wolski
2012-10-26
编辑：Sean de Wolski
2012-10-26
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
2012-10-27
@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.
更多回答（5 个）
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 个评论
Azzi Abdelmalek
2012-10-26
编辑：Azzi Abdelmalek
2012-10-26
I am not getting the result
Matt J
2012-10-26
Works fine for me.
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
2012-10-26
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
2012-10-26
编辑：Azzi Abdelmalek
2012-10-26
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.
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
2012-10-26
@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
2012-10-26
Yes, my tes is'nt adequate, like you said, I must make a test in a loop.
Azzi Abdelmalek
2012-10-26
eval(sprintf('isPeak(start_peak_hour(%d):end_peak_hour(%d)) = 1,',[1:n_peak;1:n_peak]))
3 个评论
Sam
2012-10-26
Azzi Abdelmalek
2012-10-26
If you want to avoid eval
isPeak(cell2mat(arrayfun(@(x,y) x:y, start_peak_hour,end_peak_hour,'un',0)))=1
Jan
2012-10-26
And you should avoid EVAL...
[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);
3 个评论
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
2012-10-26
At least I've considered the last element...
True!
Chris A
2012-10-26
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'),
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
类别
在 帮助中心 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
