Re: I NEED HELP: Differences between MATLAB and Octave

[Paul Kienzle]

On Fri, Jun 07, 2002 at 02:25:17PM -0500, Marco Boni wrote:
> The simple code I've attached to this mail runs very fast under Matlab,
> while it is executed very slow under Octave.
> 
> Why? Where I've failed!

For-loops in octave are embarrassingly slow.  You have nested
forloops (one in your code, and one in hist).  Both are probably
removable with a bit of cleverness.

For example, in hist:

        for i = 1:n-1
                cutoff (i) = (x (i) + x (i + 1)) / 2;
        endfor

is simply translated into

        cutoff = ( x(1:n-1) + x(2:n) ) / 2;

The next loop in hist.m is more tricky:

        freq(1) = sum(y < cutoff (1));
        for i = 2:n-1
                freq (i) = sum (y >= cutoff (i-1) & y < cutoff (i));
        endfor
        freq(n) = sum (y >= cutoff (n-1));

You can get part way there with sort:

        [s, idx] = sort ( [cutoff(:); y(:)] );
        pt = [0; idx > n; 0];

This gives you a sequence of 0's and 1's where the zeros represent
bin boundaries and the ones represent bin contents.  Because sort
is 'stable' all y values which are identical to an x value will
appear after the x value, so the semi-open [x-1, x) logic will
be preserved.  The leading and trailing zeros capture all the y's
which are not otherwise in a bin.

Using cumulative sum, you can turn the ones into frequency counts

        chist = cumsum(pt); 

But you only need the counts at the bin boundaries:

        chist = s(find(diff(chist) == 0));

Now you have a 'cumulative histogram'.  Differentiate it and you 
should have the histogram:

        freq = diff(chist);

I may messed up something along the way, but hopefully you get
the idea.  Could you fix it, test it and submit it to bug-octave?

Since you are dealing with really large numbers of bins (up to 2^20 if
I'm reading your code correctly), fixing hist should address most of
the speed problems.  

You can also remove a number of statements outside of your for loop,
but that probably won't gain you very much.

Paul Kienzle
address@hidden

> 
> 
> 
> --
> Marco Boni
> Viale delle Magnolie, 4
> 50142 - Firenze
> Cell.: +39 335 6079353
> Home : +39 055 7398033
> Email: address@hidden
> -----------------------------------------------------
>  My machine is an i686 and running Red Hat Linux 7.2
>                Why don't you try it?
> -----------------------------------------------------
> 
> 

> %entropia di un segnale random
> 
> closeplot;
> clear all;
> k=16;
> nt=2^k;
> a=floor(rand(1,nt)+.5);
> dt=k+4;
> for n=1:dt
>    dimb=n;
>    numb=floor(nt/dimb);
>    dim=dimb*numb;
>    aa=reshape(a(1:dim),dimb,numb);
>    x=(1:1:n)-1;
>    xx=2.^x;
>    ris1=xx*aa+1;
>    ris2=hist(ris1,2^n);
>    ris3=ris2(find(ris2));
>    %
> size(ris3,2)/size(ris2,2);
>    n
>    p=ris3/numb;
>    ris(n)=-sum(p.*log(p));
> end;
> entr=diff(ris);
> entr(1)=log(2);
> plot(entr);hold;plot(entr,'*');
> plot(log(2)*ones(1,dt),'r');grid
> input ("Pick a number, any number! ")



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