From: Lex on
Dear all,

how can I vectorize this pice of code? Let`s say m=100.
Thanks for your help.

[m n] = size(x);
y = sort(x);
for i=1:m
for j=1:m
ecop(i,j) = sum( (x(:,1)<=y(i,1)).*(x(:,2)<=y (j,2)) )/m;
end
end
From: Miroslav Balda on
Lex <laux5101(a)gmx.net> wrote in message <1856195661.15616.1277719543151.JavaMail.root(a)gallium.mathforum.org>...
> Dear all,
>
> how can I vectorize this pice of code? Let`s say m=100.
> Thanks for your help.
>
> [m n] = size(x);
> y = sort(x);
> for i=1:m
> for j=1:m
> ecop(i,j) = sum( (x(:,1)<=y(i,1)).*(x(:,2)<=y (j,2)) )/m;
> end
> end

Hi Lex,

I am not sure if your code can be vectorized, nevertheless, it could be accelerated by simple modifications. The following code contains a slightly modified original, that improves behavior of the code by pre-definition of ecop matrix, and its improvement in the second half:

% Lex.m 2010-06-28
%%%%%%%%
m = inp('m',500,'%5d'); % FEX Id: 9033
x = rand(m,2);
y = sort(x);
ecop = zeros(m);
e = zeros(m);


tic
for i=1:m
for j=1:m
ecop(i,j) = sum( (x(:,1)<=y(i,1)).*(x(:,2)<=y (j,2)) );
end
end
ecop = ecop/m;
toc

tic
for i=1:m
row = single(x(:,1)<=y(i,1))';
for j=1:m
e(i,j) = row * single(x(:,2)<=y(j,2));
end
end
e = e/m;
toc

normdel = norm(ecop-e)

The computing times of the second part are about one half of the first ones.
Best regards,

Mira
From: Matt Fig on
m = length(x);
y = sort(x);
AA = bsxfun(@le,x(:,1),y(:,1)');
BB = bsxfun(@le,x(:,2),y(:,2)');
ecop = double(AA')*double(BB)/m;
From: Bruno Luong on
Here is a fast one (using the following file from FEX)
http://www.mathworks.com/matlabcentral/fileexchange/23897-n-dimensional-histogram

% Data
x=ceil(50*rand(1000,2));

% OP's for-loop engine
tic
[m n] = size(x);
y = sort(x);
ecop1 = zeros(m);
for i=1:m
for j=1:m
ecop(i,j) = sum( (x(:,1)<=y(i,1)).*(x(:,2)<=y (j,2)) )/m;
end
end
toc % 29.341238 seconds.

% Matt's engine
tic
[m n] = size(x);
y = sort(x);
AA = bsxfun(@le,x(:,1),y(:,1)');
BB = bsxfun(@le,x(:,2),y(:,2)');
ecop2 = double(AA')*double(BB)/m;
toc % 0.309713 seconds.

% Bruno's HISTCN engine
tic
[m n] = size(x);
y = sort(x);
[e1 I1 J1]=unique(y(:,1),'last');
[e2 I2 J2]=unique(y(:,2),'last');
e = cumsum(cumsum(histcn(x, e1, e2),1),2);
ecop3 = e(J1,J2)/m;
toc % Elapsed time is 0.036968 seconds.

% Bruno
From: Lex on
Dear all,

thank you very much for your suggestions. This helped me a lot.

Kind regards, Pat