Re: How do I replace this "for" loop?

[David Bateman]

Tom Holroyd wrote:

> Oops!  There was a(nother) bug in my script -- the original line was
>
>         x = 2*(x>p);
>
> but it should have been
>
>         x = 2*x(x>p);
>
> and the really astonishing part is, the second line is faster.
> In fact (see attached graph), it shows that in this case, find is 
> *never* faster than using logical indexing.  Not only that, even the 
> logical indexing gets faster as the number of elements selected 
> decreases!  I'm not sure why that is.
There are some further provisos, that are due to certain inefficiencies 
in the octave indexing code as it stands... When you do the "x>p" 
operation you create a boolean array, where a boolean in C++ is a one 
byte element, so if x has n elements you use n bytes. Furthermore, you 
then have to consider how the idx_vector class treats logical indexing. 
The internal representation of idx_vector is of the type octave_idx_type 
which in general is a 4-byte integer (though can be 8-bytes in some 
circumstances with recent octave versions). Look at the constructor of 
an idx_vector from a boolean array.

IDX_VEC_REP::idx_vector_rep (const boolNDArray& bnda)
 : data (0), len (bnda.length ()), num_zeros (0), num_ones (0),
   max_val (0), min_val (0), count (1), frozen_at_z_len (0),
   frozen_len (0), colon (0), one_zero (1), initialized (0),
   frozen (0), colon_equiv_checked (0), colon_equiv (0),
   orig_dims (bnda.dims ())
{
 if (len == 0)
   {
     initialized = 1;
     return;
   }
 else
   {
     octave_idx_type k = 0;
     data = new octave_idx_type [len];

     for (octave_idx_type i = 0; i < len; i++)
   data[k++] = tree_to_mat_idx (bnda.elem (i));
   }

 init_state ();
}

so you create a direct copy of the boolean object in the type 
octave_idx_type and therefore use either 4*n or 8*n bytes. The 
assign/index functions then explicitly check the flag one_zero (above 
set to 1) to see if the indexing is logical or not... The problem is 
that the indexing therefore takes comparable memory to the original 
matrix, and in some circumstances much more (eg indexing a logical or 
int8 array). So if you are working at the limit of your available memory 
this is bad.

I see there are two ways to address this issue

1) Make the above constructor handle the conversion from logical 
indexing to a normal indexing which if there are fewer than 25% of the 
elements of the logical indexing at true will take only the same amount 
of memory as the indexing itself. But there is still a cost in terms of
memory use. The big advantage is that there is no code changes needed in 
the assign/indexing functions, and a further advantage is that
we could add the constructor for SparseBoolMatrix and something like

a = sprandn (1e6, 1e6, 1e5);
idx = a < 0.;
a(idx) = 0;

then won't explode the memory as it currently does (it does in matlab as 
well).

2) Allow the idx_vector class store the logical array. This will be 
cheaper in terms of the memory use, but the () operator would need to be 
adjusted to select where the data comes from with a potential slow up of 
all indexing operators. So this option is probably not the right one. 
Also the sparse boolean indexing case above would mean that there would 
have to be multiple means of storing the indexing..

Sorry, I started a thread about the internals of indexing, but after 
this I have the impression that 1) above won't be that difficult and the 
pay-off in terms of memory and perhaps even speed might make this a 
worthwhile change to try..

Cheers
David

--
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax) 
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary



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