Fwd: Re: Re: Another vectorisation question

[John Reed]

>>>> Jaroslav Hajek <address@hidden> 7/15/2010 11:00 PM >>>
> On Fri, Jul 16, 2010 at 3:24 AM, John Reed <address@hidden> wrote:
>> Hello,
>>
>> I am looking at a nearest neighbor problem.  I have a rectangular grid and 
>> random coordinates that fall within the 2d grid.  I want to create an array 
>> that stores all the random coordinates that are nearest to each respective 
>> grid point.
>>
>> I've taken two different approaches to solving this problem, one uses the 
>> find function, and the other the histc function, but in either case I don't 
>> know how to get rid of all the for loops and i would like to speed up the 
>> code.
>>
>> Is there a way to vectorise either of these approaches?  or is there some 
>> totally different approach that would give me the answer I want without for 
>> loops?
>>
>> The two different approaches I've tried are below.
>>
>> Any help would be greatly appreciated,
>>
>> John
>> ..........
>>
>> incr = 70;   % maximum x & y integer value on square unit grid
>> pt_num = 40000;
>> coord=incr*rand(pt_num,2);
>>
>> % 1st Approach - find grid point that coordinate is closest to and assign to 
>> array
>> % holding all coordinates closest to each respective grid point
>> tic;
>> for i = 0:incr
>>   for j = 0:incr
>>      x = [ i - 0.5; i + 0.5];
>>      y = [ j - 0.5; j + 0.5];
>>
>>      indices = find( coord(:,1) >= x(1) & ...
>>                   coord(:,1) <  x(2) & ...
>>                   coord(:,2) >= y(1) & ...
>>                   coord(:,2) <  y(2) );
>>
>> % build array storing random coordinates indices associated with each grid 
>> coordinate
>>      num_pt_per_grid(i+1,j+1)=length(indices);
>>      if ( length(indices) >= 1 )
>>         coords_to_grid(i+1,j+1,1:length(indices)) = indices;
>>      end
>>   end
>> end
>> toc
>>
>> fflush(stdout);
>>
>> % 2nd approach - using histc group all coordinates closest to grid point and 
>> assign
>> % to array holding all coordinates closest to each grid point
>> tic;
>> % create edges for histogram
>> histo_edges = [-0.5:1:incr+.5];
>>
>> [x_count, x_index] = histc(coord(:,1),histo_edges);
>> [y_count, y_index] = histc(coord(:,2),histo_edges);
>>
>> % build array storing random coordinates indices associated with each grid 
>> coordinate
>> num_pt_per_grid2 = zeros(incr+1,incr+1);
>> for i=1:length(x_index)
>>   num_pt_per_grid2(x_index(i),y_index(i))++;
>>   
>> coords_to_grid2(x_index(i),y_index(i),num_pt_per_grid2(x_index(i),y_index(i)))
>>  = i;
>> end
>> toc
>>
>>
>>
>
>
> % 3rd approach - use a cell array of vectors to accumulate.
> tic;
>
> edges = [-0.5:1:incr+.5];
> idx = lookup (edges, coord);
>
> coords_to_grid3 = accumarray (idx, 1:pt_num, [incr, incr]+1, @(x) {x});
> toc
>
>
> now, the indices falling into grid (i,j) are given by
> coords_to_grid3{i,j} instead of coords_to_grid2(i,j,1:num_pt_per_grid(i,j)).
>
>
> timings in 3.3.51+ (recommended for best performance):
>
> octave:1> tt
> Elapsed time is 2.06051 seconds.
> Elapsed time is 0.955433 seconds.
> Elapsed time is 0.00869894 seconds.
>
> coords_to_grid3 can be converted to coords_to_grid2 with a little
> extra effort (if that is needed at all).
>
> hth
>
> --
> RNDr. Jaroslav Hajek, PhD
> computing expert & GNU Octave developer
> Aeronautical Research and Test Institute (VZLU)
> Prague, Czech Republic
> url: www.highegg.matfyz.cz 
>
>

Thank you very much for the suggestion Jaroslav.

I'll test it out in my code and try to learn more about cell arrays.

John