Re: de-for-ing..

[Miroslaw Kwasniak]

On Sun, Feb 06, 2005 at 10:42:36PM +0100, Gorazd Brumen wrote:
> The results are the same, but the for loop is much
> faster and less memory consuming. the kron solution
> is preety much stupid.

I'm suprized too, because octave has slow loops.

> 
> Is there a better way than the for loop?

I've tested also two additional methods with repmat: and sort:

place3=sum(repmat(x_init(:),1,lx) <= repmat(x(:)',n,1));
t3=cputime;

and sort:

[y i0]=sort([x_init, x]);
i=find(i0>n);
place4(i0(i)-n)=cumsum([i(1) diff(i)]-1);
t4=cputime;

Results for
ftest(n=1000, lx=10000);
ftest(n=10000, lx=1000);

cpu-tme in seconds (2.1.63 on A1600XP)
  for-loop   kronecker  repmat     sort
  1.600000   2.350000   4.200000   0.040000
  0.8700000  2.3800000  2.5200000  0.0100000

Winner is the method with sort, but it may fail, because it assumes that
sort doesn't change order for equal values:

   x(i) == x(j) & i < j ==> find(i0==i) < find(i0==j)

and it may depend on sort implementation.


Mirek

Test procedure:
===============================
0;

function [place1,place4,i,x_init,x]=ftest(n,lx)
if length(n)>1
  x_init=n; n=length(n);
  x=lx; lx=length(lx);
else
  x_init=(1:n)/n*lx;
  x=1:lx;
endif

t0=cputime;
place1=zeros(1,lx);
for i = 1:lx
    place1(i) = sum( x_init <= x(i) );
end
t1=cputime;

x_ext = kron (x, ones (1,n));
x_init_ext = kron (ones (1,lx), x_init );
x_tog = (x_init_ext <= x_ext);
place2=sum(reshape (x_tog,n,lx));
t2=cputime;

place3=sum(repmat(x_init(:),1,lx) <= repmat(x(:)',n,1));
t3=cputime;

[y i0]=sort([x_init, x]);
i=find(i0>n);
place4(i0(i)-n)=cumsum([i(1) diff(i)]-1);
t4=cputime;

disp(diff([t0 t1 t2 t3 t4]));
disp([any(place1~=place2) any(place1~=place3)  any(place1~=place4)]);

endfunction

n=100
ftest(10*n, 100*n);
ftest(100*n, 10*n);
ftest(randn(1,10*n), randn(1,100*n));
[q1,q4,qi,xi,x]=ftest(randn(1,100*n), randn(1,10*n));



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------