清空筛选条件
清空筛选条件

vectorize a for loop

1 次查看（过去 30 天）
Sam
Sam 2012-10-26
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
José-Luis 2012-10-26
Loops are not always evil, and in this case probably even faster than the alternatives.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Sean de Wolski
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 个评论
显示 4更早的评论
Sean de Wolski
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
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 个）

Matt J
Matt J 2012-10-26
编辑：Matt J 2012-10-26
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 个评论
显示 6更早的评论
Sean de Wolski
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
Azzi Abdelmalek 2012-10-26
Yes, my tes is'nt adequate, like you said, I must make a test in a loop.

请先登录，再进行评论。

Azzi Abdelmalek
Azzi Abdelmalek 2012-10-26
eval(sprintf('isPeak(start_peak_hour(%d):end_peak_hour(%d)) = 1,',[1:n_peak;1:n_peak]))
  3 个评论
显示 1更早的评论
Azzi Abdelmalek
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
Jan 2012-10-26
And you should avoid EVAL...

请先登录，再进行评论。

Jan
Jan 2012-10-26
编辑：Jan 2012-10-27
[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 个评论
显示 1更早的评论
Jan
Jan 2012-10-26
编辑：Jan 2012-10-26
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
Matt J 2012-10-26
At least I've considered the last element...
True!

请先登录，再进行评论。

Chris A
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'),
Jan
Jan 2012-10-27
编辑：Jan 2012-10-27
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

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by