From help-request at octave dot org Mon Feb  7 07:15:53 2005
Subject: Re: de-for-ing..
From: Miroslaw Kwasniak <mirek at zind dot ikem dot pwr dot wroc dot pl>
To: help at octave dot org
Date: Mon, 7 Feb 2005 14:20:29 +0100

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
-------------------------------------------------------------