Thread Subject: loop & find are slow

dear all,

i have the following piece of code in my function:

for iipixel=1:numberpixels
    % do something...
    pixels_idx=find(row_idx_matrix==row(iipixel) &
col_idx_matrix==scol(iipixel));
    % continue
end

when i put on the profiler on i see that these lines are
making the code really slow (90% of the time is spend
there)

i have to find the value that fullfil 2 conditions.
row_idx_matrix and col_idx_matrix are 2 matrices with a
size of about 250 by 250

is there any way to make this go faster?

thanks,
carmelo

"carmelo " <km.nospam.scriba.nospam@yahoo.esnospam> writes:

> dear all,
>
> i have the following piece of code in my function:
>
> for iipixel=1:numberpixels
> % do something...
> pixels_idx=find(row_idx_matrix==row(iipixel) &
> col_idx_matrix==scol(iipixel));
> % continue
> end
>
> when i put on the profiler on i see that these lines are
> making the code really slow (90% of the time is spend
> there)
>
> i have to find the value that fullfil 2 conditions.
> row_idx_matrix and col_idx_matrix are 2 matrices with a
> size of about 250 by 250

A histogram function might work here, but it is hard to say without more
information.

Are the matrices integer-valued? Over what range? What value is
numberpixels? What do you do with the resulting answers? There might
be a way to vectorize the whole loop.

-Peter

Peter Boettcher <boettcher@ll.mit.edu> wrote in message
<muyve8vp4ab.fsf@G99-Boettcher.llan.ll.mit.edu>...
> "carmelo " <km.nospam.scriba.nospam@yahoo.esnospam>
writes:
>
> > dear all,
> >
> > i have the following piece of code in my function:
> >
> > for iipixel=1:numberpixels
> > % do something...
> > pixels_idx=find(row_idx_matrix==row(iipixel) &
> > col_idx_matrix==scol(iipixel));
> > % continue
> > end
> >
> > when i put on the profiler on i see that these lines
are
> > making the code really slow (90% of the time is spend
> > there)
> >
> > i have to find the value that fullfil 2 conditions.
> > row_idx_matrix and col_idx_matrix are 2 matrices with
a
> > size of about 250 by 250
>
> A histogram function might work here, but it is hard to
say without more
> information.
>
> Are the matrices integer-valued?
Yes, they are

> Over what range?
1 to 300

> What value is numberpixels?
numberpixels is about 40000

> What do you do with the resulting answers?
I use pixels_idx as an index to select a few values from
another matrix

> There might be a way to vectorize the whole loop.
> -Peter

Please let me know if you need anything else.
Thanks in advance!!

"carmelo " <km.nospam.scriba.nospam@yahoo.esnospam> writes:

> Peter Boettcher <boettcher@ll.mit.edu> wrote in message
> <muyve8vp4ab.fsf@G99-Boettcher.llan.ll.mit.edu>...
>> "carmelo " <km.nospam.scriba.nospam@yahoo.esnospam>
> writes:
>>
>> > dear all,
>> >
>> > i have the following piece of code in my function:
>> >
>> > for iipixel=1:numberpixels
>> > % do something...
>> > pixels_idx=find(row_idx_matrix==row(iipixel) &
>> > col_idx_matrix==scol(iipixel));
>> > % continue
>> > end
>> >
>> > when i put on the profiler on i see that these lines
> are
>> > making the code really slow (90% of the time is spend
>> > there)
>> >
>> > i have to find the value that fullfil 2 conditions.
>> > row_idx_matrix and col_idx_matrix are 2 matrices with
> a
>> > size of about 250 by 250
>>
>> A histogram function might work here, but it is hard to
> say without more
>> information.
>>
>> Are the matrices integer-valued?
> Yes, they are
>
>> Over what range?
> 1 to 300

So row, row_idx_matrix, scol, and col_idx_matrix all have values that
range from 1 to 300?

>> What value is numberpixels?
> numberpixels is about 40000

Are the row/col pairs unique?

>> What do you do with the resulting answers?
> I use pixels_idx as an index to select a few values from
> another matrix

And what else...

>> There might be a way to vectorize the whole loop.
>> -Peter
>
> Please let me know if you need anything else.

Lots more! What is your loop doing? Start at the top level, and
describe the problem. If you can, boil down your loop code to the bare
minimum and actually post it.

-Peter

Peter Boettcher <boettcher@ll.mit.edu> wrote in message
<muyhckfowce.fsf@G99-Boettcher.llan.ll.mit.edu>...
> "carmelo " <km.nospam.scriba.nospam@yahoo.esnospam>
writes:
>
> > Peter Boettcher <boettcher@ll.mit.edu> wrote in
message
> > <muyve8vp4ab.fsf@G99-Boettcher.llan.ll.mit.edu>...
> >> "carmelo " <km.nospam.scriba.nospam@yahoo.esnospam>
> > writes:
> >>
> >> > dear all,
> >> >
> >> > i have the following piece of code in my function:
> >> >
> >> > for iipixel=1:numberpixels
> >> > % do something...
> >> > pixels_idx=find(row_idx_matrix==row(iipixel) &
> >> > col_idx_matrix==scol(iipixel));
> >> > % continue
> >> > end
> >> >
> >> > when i put on the profiler on i see that these
lines
> > are
> >> > making the code really slow (90% of the time is
spend
> >> > there)
> >> >
> >> > i have to find the value that fullfil 2 conditions.
> >> > row_idx_matrix and col_idx_matrix are 2 matrices
with
> > a
> >> > size of about 250 by 250
> >>
> >> A histogram function might work here, but it is hard
to
> > say without more
> >> information.
> >>
> >> Are the matrices integer-valued?
> > Yes, they are
> >
> >> Over what range?
> > 1 to 300
>
> So row, row_idx_matrix, scol, and col_idx_matrix all
have values that
> range from 1 to 300?
>
> >> What value is numberpixels?
> > numberpixels is about 40000
>
> Are the row/col pairs unique?
>
> >> What do you do with the resulting answers?
> > I use pixels_idx as an index to select a few values
from
> > another matrix
>
> And what else...
>
> >> There might be a way to vectorize the whole loop.
> >> -Peter
> >
> > Please let me know if you need anything else.
>
> Lots more! What is your loop doing? Start at the top
level, and
> describe the problem. If you can, boil down your loop
code to the bare
> minimum and actually post it.
>
> -Peter


I am sorry if i did not provide enought details. i was
just trying to simplify the problem and the loop.
here you have more details

i have 2 matrices. one represents a classified high
spatial resolution image and another one a low spatial
resolution image. they images are over the same area so
one low spatial resolution pixel has several classified
high spatial resolution pixels "inside of it".

what i try to do with this loop is desaggregate the low
spatial resolution pixel values:

% initialize output1 and output2 and set upper and lower
boundaries for the solution of the lsqlin problem

for iipixel=1:numberpixels
    % prepare PV (m x1) and F(mxn) matrices.
    % iipixel is used as a pointer to extract PV and F
from 2 different
    % matrices
    PV=extract_pixel_values
(MatrixLOWpadded,windowsize,paddedrow(iipixel),paddedcol
(iipixel));
    [F classnumber]=extract_fractions
(Fractionmatrix,windowsize,paddedrow(iipixel),paddedcol
(iipixel));
    
    % solving a lsqlin problem with constraints
    endmembers=sls_alloc(F,PV,lb,ub); %sls_alloc is not a
standard function in matlab
    RMSE=sqrt(sum(((PV-F*endmembers).^2))/length(PV)); %
compute Root Mean Square error

    % Finding the indices of the high spatial resolution
pixels under each low spatial resolution pixel
    pixels_idx=find(row_idx_matrix==row(iipixel) &
col_idx_matrix==col(iipixel));
    
    % Replace classnumbers by their corresponding values
(endmembers)
    classes_inLOWPIXEL=classified_image(pixels_idx);
    endmembers_inLOWPIXEL=classesinLOWPIXEL*0;
    for jj=1:length(classnumber)
        endmembers_inLOWPIXEL
(classes_inLOWPIXEL==classnumber(jj))=endmembers(jj);
    end
    
    % Store outputs
    output1(pixels_idx)=endmembers_inLOWPIXEL;
    output2(row(iipixel),col(iipixel))=RMSE;
end


row_idx_matrix and col_idx_matrix are 2 matrices with the
size of the high spatial resolution matrix that contain
the row and col of their corresponding low spatial
resolution pixel. Because the low spatial resolution image
is about 300 by 300 pixels these matrices have values
between 1 and 300 (aprox)

> Are the row/col pairs unique?
Yes, they are all unique. i use iipixel as an index
because the low spatial resolution image has "mask" (zero
values) and like this i select the row/col pairs with
nonzero values.

Please let me know if this is clear (i am afraid that i
was not very clear... it is difficult for me to boil down
the loop and to explain the problem to an "outsider" )
hopefully, it is a bit clearer now.

many thanks
Carmelo

"carmelo " <km.nospam.scriba.nospam@yahoo.esnospam> writes:
> I am sorry if i did not provide enought details. i was
> just trying to simplify the problem and the loop.
> here you have more details

Simplifying is good. Thanks for trying! I was hoping there was a way
to vectorize the whole loop, but the whole thing appears very complex,
so I think you probably got the right level of simplifying after all.

> i have 2 matrices. one represents a classified high
> spatial resolution image and another one a low spatial
> resolution image. they images are over the same area so
> one low spatial resolution pixel has several classified
> high spatial resolution pixels "inside of it".
>
> what i try to do with this loop is desaggregate the low
> spatial resolution pixel values:
>
> % initialize output1 and output2 and set upper and lower
> boundaries for the solution of the lsqlin problem
>
> for iipixel=1:numberpixels
> % prepare PV (m x1) and F(mxn) matrices.
> % iipixel is used as a pointer to extract PV and F
> from 2 different
> % matrices
> PV=extract_pixel_values
> (MatrixLOWpadded,windowsize,paddedrow(iipixel),paddedcol
> (iipixel));
> [F classnumber]=extract_fractions
> (Fractionmatrix,windowsize,paddedrow(iipixel),paddedcol
> (iipixel));
>
> % solving a lsqlin problem with constraints
> endmembers=sls_alloc(F,PV,lb,ub); %sls_alloc is not a
> standard function in matlab
> RMSE=sqrt(sum(((PV-F*endmembers).^2))/length(PV)); %
> compute Root Mean Square error
>
> % Finding the indices of the high spatial resolution
> pixels under each low spatial resolution pixel
> pixels_idx=find(row_idx_matrix==row(iipixel) &
> col_idx_matrix==col(iipixel));
>
> % Replace classnumbers by their corresponding values
> (endmembers)
> classes_inLOWPIXEL=classified_image(pixels_idx);
> endmembers_inLOWPIXEL=classesinLOWPIXEL*0;
> for jj=1:length(classnumber)
> endmembers_inLOWPIXEL
> (classes_inLOWPIXEL==classnumber(jj))=endmembers(jj);
> end
>
> % Store outputs
> output1(pixels_idx)=endmembers_inLOWPIXEL;
> output2(row(iipixel),col(iipixel))=RMSE;
> end
>
>
> row_idx_matrix and col_idx_matrix are 2 matrices with the
> size of the high spatial resolution matrix that contain
> the row and col of their corresponding low spatial
> resolution pixel. Because the low spatial resolution image
> is about 300 by 300 pixels these matrices have values
> between 1 and 300 (aprox)

Sorry, I still can't follow the indexing used. When you talk about
high-resolution pixels "corresponding" to low-res pixels, I think that
the relationship might be regular. That is, pixel 1 low-res
corresponds to pixels 1 to 300 (or whatever) high-res. If that's so,
you don't need a "find" at all.

I will try to explain my suggestion in general terms:

Twice in your code you take a matrix with a relatively few number of
integer values. You loop through all the values that the matrix
contains, find the indices of that value, and set the result to
something. The first time is big pixels_idx find. The second one is
the endmembers_inLOWPIXEL. Try to convert these instead to a look-up
table, where each entry in the matrix is examined ONCE, and the
appropriate value is stored.

Take this example (I have no idea if the sizes/shapes are right, but you
get the idea):

-------------------------------
classnumber = [1:100];
endmembers = rand(1,100);
classes_inLOWPIXEL = ceil(100*rand(1,100000));

tic
for jj=1:length(classnumber)
  endmembers_inLOWPIXEL(classes_inLOWPIXEL==classnumber(jj))=endmembers(jj);
end
toc

tic
  endmembers_inLOWPIXEL2 = endmembers(classes_inLOWPIXEL);
toc

isequal(endmembers_inLOWPIXEL, endmembers_inLOWPIXEL2)
---------------------------------

Elapsed time is 0.076110 seconds.
Elapsed time is 0.002842 seconds.


Do you see the difference? Instead of looping over all values, use them
as a lookup. I realize this particular part of your code is not a
bottleneck. But I think maybe the same concept could be applied to
building output1. For instance, classes_inLOWPIXEL is just a subset of
classified_image, and the output only goes into output1. So already we
can do:

output1(pixels_idx) = endmembers(classified_image(pixels_idx));

The part you'll have to work on is now the outer loop. Perhaps you can
build a new LUT by modifying the values of classified_image. Think
about some code which leads to this structure:

output1 = endmembers_lut(classified_image_modified);

The trick is to somehow NOT look at each element of row_idx_matrix and
col_idx_matrix many times.


-Peter