kdtree2.lyx

kd tree java implementation 2

\newline 
    ! Find the 'nn' vectors in the tree nearest to 'qv' in euclidean norm
\newline 
    ! returning their indexes and distances in 'indexes' and 'distances'
\newline 
    ! arrays already allocated passed to this subroutine.
\newline 
    type (kdtree2), pointer      :: tp
\newline 
    real, target, intent (In)    :: qv(:)
\newline 
    integer, intent (In)         :: nn
\newline 
    type(kdtree2_result), target :: results(:)
\newline 

\newline 
subroutine kdtree2_n_nearest_around_point(tp,idxin,correltime,nn,results)
\newline 
    ! Find the 'nn' vectors in the tree nearest to point 'idxin',
\newline 
    ! with correlation window 'correltime', returing results in
\newline 
    ! results(:), which must be pre-allocated upon entry.
\newline 
    type (kdtree2), pointer        :: tp
\newline 
    integer, intent (In)           :: idxin, correltime, nn
\newline 
    type(kdtree2_result), target   :: results(:)
\newline 

\newline 
subroutine kdtree2_r_nearest_around_point(tp,idxin,correltime,r2,nfound,nalloc,results)
\newline 
    type (kdtree2), pointer      :: tp
\newline 
    integer, intent (In)         :: idxin, correltime, nalloc
\newline 
    real, intent(in)             :: r2
\newline 
    integer, intent(out)         :: nfound
\newline 
    type(kdtree2_result), target :: results(:)
\newline 

\newline 
function kdtree2_r_count(tp,qv,r2) result(nfound)
\newline 
    ! Count the number of neighbors within square distance 'r2'. 
\newline 
    type (kdtree2), pointer   :: tp
\newline 
    real, target, intent (In) :: qv(:)
\newline 
    real, intent(in)          :: r2
\newline 
    integer                   :: nfound
\newline 

\newline 
function kdtree2_r_count_around_point(tp,idxin,correltime,r2)  result(nfound)
\newline 
    type (kdtree2), pointer :: tp
\newline 
    integer, intent (In)    :: correltime, idxin
\newline 
    real, intent(in)        :: r2
\newline 
    integer                 :: nfound
\newline 

\newline 
subroutine kdtree2_n_nearest_brute_force(tp,qv,nn,results) 
\newline 
    ! find the 'n' nearest neighbors to 'qv' by exhaustive search.
\newline 
    ! only use this subroutine for testing, as it is SLOW!  The
\newline 
    ! whole point of a k-d tree is to avoid doing what this subroutine
\newline 
    ! does.
\newline 
    type (kdtree2), pointer :: tp
\newline 
    real, intent (In)       :: qv(:)
\newline 
    integer, intent (In)    :: nn
\newline 
    type(kdtree2_result)    :: results(:) 
\newline 

\newline 
subroutine kdtree2_sort_results(nfound,results)
\newline 
    !  Use after search to sort results(1:nfound) in order of increasing 
\newline 
    !  distance.
\newline 
    integer, intent(in)          :: nfound
\newline 
    type(kdtree2_result), target :: results(:) 
\newline 

\backslash 
end{verbatim}
\end_inset 


\layout Standard
\paragraph_spacing single 
\noindent 
An example follows.
\begin_inset ERT
status Collapsed

\layout Standard

\backslash 
begin{verbatim}
\newline 
program kdtree2_example
\newline 
  use kdtree2_module
\newline 
  type(kdtree2), pointer            :: tree
\newline 
  integer                           :: N,d
\newline 
  real, allocatable                 :: mydata(:,:)
\newline 
  type(kdtree2_result), allocatable :: results(:)
\newline 

\newline 
! user sets d, the dimensionality of the Euclidean space
\newline 
! and N, the number of points in the set. 
\newline 

\newline 
allocate(mydata(d,N))   
\newline 
! note order, d is first, N second. 
\newline 

\newline 
! read in vectors into mydata(j,i) for j=1..d, i=1..N
\newline 

\newline 
tree => kdtree2_create(mydata,rearrange=.true.,sort=.true.)
\newline 
! Create the tree, ask for internally rearranged data for speed,
\newline 
! and for output sorted by increasing distance from the
\newline 
! query vector
\newline 

\newline 
allocate(results(20))
\newline 
call kdtree2_n_nearest_around_point(tree,idxin=100,nn=20,correltime=50,results)
\newline 

\newline 
! Now the 20 nearest neighbors to mydata(*,100) are in results(:) except
\newline 
! that points within 50 time units of idxin=50 are not considered as valid neighbors.
\newline 
!
\newline 
write (*,*) 'The first 10 near neighbor distances are: ', results(1:10)%dis
\newline 
write (*,*) 'The first 10 near neighbor indexes   are: ', results(1:10)%idx
\newline 

\backslash 
end{verbatim}
\end_inset 


\layout Subsection

C++
\layout Standard
\paragraph_spacing single 
\noindent 
The interface header is 
\family typewriter 
kdtree2.hpp
\family default 
 and main code in 
\family typewriter 
kdtree2.cpp
\family default 
.
 The BOOST (
\family typewriter 
www.boost.org
\family default 
) library must be installed
\begin_inset Foot
collapsed false

\layout Standard

On the author's Fedora Core 2 Linux system, this can be done by installing
 the 
\family typewriter 
boost
\family default 
 and 
\family typewriter 
boost-devel
\family default 
 RPM packages.
\end_inset 

 as should the Standard Template library.
 Interfaces for important public routines follow.
 Note that sorting of results in increasing distance can by done using STL
 as 
\family typewriter 
sort(results.begin(),results.end())
\family default 
.
\layout Standard
\paragraph_spacing single 
\noindent 

\begin_inset ERT
status Collapsed

\layout Standard

\backslash 
begin{verbatim}
\newline 
  //constructor
\newline 
  kdtree2(kdtree2_array& data_in,bool rearrange_in = true,int dim_in=-1);
\newline 
  // destructor
\newline 
  ~kdtree2();
\newline 
  // set to true to always sort
\newline 
  bool sort_results;
\newline 

\newline 
  void n_nearest(vector<float>& qv, int nn, kdtree2_result_vector& result);
\newline 
  // search for n nearest to a given query vector 'qv'.
\newline 
  void n_nearest_around_point(int idxin, int correltime, int nn,
\newline 
      kdtree2_result_vector& result);
\newline 
  // search for 'nn' nearest to point [idxin] of the input data, excluding
\newline 
  // neighbors within correltime 
\newline 
  void r_nearest(vector<float>& qv, float r2,kdtree2_result_vector& result); 
\newline 
  // search for all neighbors in ball of size (square Euclidean distance)
\newline 
  // r2.   Return number of neighbors in 'result.size()', 
\newline 
  void r_nearest_around_point(int idxin, int correltime, float r2,
\newline 
      kdtree2_result_vector& result);
\newline 
  // like 'r_nearest', but around existing point, with decorrelation
\newline 
  // interval. 
\newline 
  int r_count(vector<float>& qv, float r2);
\newline 
  // count number of neighbors within square distance r2.
\newline 
  int r_count_around_point(int idxin, int correltime, float r2);
\newline 
  // like r_count, but around an extant point.
\newline 

\backslash 
end{verbatim}
\end_inset 


\layout Standard
\paragraph_spacing single 
\noindent 
An example:
\begin_inset ERT
status Collapsed

\layout Standard

\backslash 
begin{verbatim}
\newline 

\newline 
#include <vector>
\newline 
#include <boost/multi_array.hpp>
\newline 

\newline 
using namespace boost;   
\newline 
using namespace std;   
\newline 

\newline 
#include "kdtree2.hpp"
\newline 

\newline 
typedef multi_array<float,2> array2dfloat; 
\newline 

\newline 
main() {
\newline 
  kdtree2               *tree;
\newline 
  int                   N,d
\newline 
  array2dfloat          mydata;
\newline 
  kdtree2_result_vector results;
\newline 

\newline 
  // user sets d, dimensionality of Euclidean space and
\newline 
  // N, number of poitns in the set.
\newline 
  
\newline 
   mydata.resize(extents[N][dim]);   
\newline 
   // get space for a N x dim matrix. 
\newline 

\newline 
   // read in vectors into mydata[i][j] for i=0..N-1, and j=0..d-1
\newline 
   // NOTE:  array is in opposite order from Fortran, and is 0-based
\newline 
   // not 1-based.   This is natural for C++ just as the other was
\newline 
   // natural for Fortran. In both cases, vectors are laid out
\newline 
   //  contiguously in memory.
\newline 

\newline 
   // notice, no need to allocate size of results, as that will be
\newline 
   // handled automatically by the STL.  results has most properties
\newline 
   // of vector<kdtree2_result>.
\newline 
  
\newline 
   tree = new kdtree2(mydata,true); // create the tree, ask to rearrange
\newline 
   tree->sort_results = true;       // sort all results.
\newline 

\newline 
   tree->n_nearest_around_point(100,50,20,results);
\newline 
   // ask for 20 nearest neighbors around point 100, with correlation window
\newline 
   // 50, push onto 'results'.
\newline 
}
\newline 

\newline 

\backslash 
end{verbatim}
\end_inset 


\layout Section

performance
\layout Standard

We now compare the performance, in searches/s, between KDTREE2 and the author's
 previous version in Fortran.
\layout Standard

First, a database of 10,000 points chosen randomly and uniformly in the
 3-d unit hypercube (in main CPU cache) query vector chosen likewise, searching
 for nearest 
\begin_inset Formula $m$
\end_inset 

 neighbors.
\layout Standard
\align center 

\begin_inset  Tabular
<lyxtabular version="3" rows="6" columns="3">
<features>
<column alignment="center" valignment="top" leftline="true" rightline="true" width="0">
<column alignment="right" valignment="top" leftline="true" width="0">
<column alignment="right" valignment="top" leftline="true" rightline="true" width="0">
<row topline="true" bottomline="true">
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\layout Standard


\begin_inset Formula $m$
\end_inset 


\end_inset 
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\layout Standard

KDTREE2
\end_inset 
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\layout Standard

old KDTREE
\end_inset 
</cell>
</row>
<row topline="true">
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\layout Standard

1
\end_inset 
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\layout Standard

415843
\end_inset 
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\layout Standard

367530
\end_inset 
</cell>
</row>
<row topline="true">
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\layout Standard

5
\end_inset 
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\layout Standard

190531
\end_inset 
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\layout Standard

160281
\end_inset 
</cell>
</row>
<row topline="true">
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\layout Standard

10
\end_inset 
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\layout Standard

127779
\end_inset 
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\layout Standard

99919
\end_inset 
</cell>
</row>
<row topline="true">
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\layout Standard

25
\end_inset 
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\layout Standard

65359
\end_inset 
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\layout Standard

41485
\end_inset 
</cell>
</row>
<row topline="true" bottomline="true">
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\layout Standard

500
\end_inset 
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\layout Standard

4794
\end_inset 
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\layout Standard

350
\end_inset 
</cell>
</row>
</lyxtabular>

\end_inset 

 
\layout Standard

For 200,000 points in 3-d unit hypercube (larger than CPU cache).
\layout Standard
\align center 

\begin_inset  Tabular
<lyxtabular version="3" rows="6" columns="3">
<features>
<column alignment="center" valignment="top" leftline="true" rightline="true" width="0">
<column alignment="right" valignment="top" leftline="true" width="0">
<column alignment="right" valignment="top" leftline="true" rightline="true" width="0">
<row topline="true" bottomline="true">
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\layout Standard


\begin_inset Formula $m$
\end_inset