PLearn 0.1
BallTreeNearestNeighbors.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // BallTreeNearestNeighbors.cc
00004 //
00005 // Copyright (C) 2004 Pascal Lamblin & Marius Muja
00006 // 
00007 // Redistribution and use in source and binary forms, with or without
00008 // modification, are permitted provided that the following conditions are met:
00009 // 
00010 //  1. Redistributions of source code must retain the above copyright
00011 //     notice, this list of conditions and the following disclaimer.
00012 // 
00013 //  2. Redistributions in binary form must reproduce the above copyright
00014 //     notice, this list of conditions and the following disclaimer in the
00015 //     documentation and/or other materials provided with the distribution.
00016 // 
00017 //  3. The name of the authors may not be used to endorse or promote
00018 //     products derived from this software without specific prior written
00019 //     permission.
00020 // 
00021 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00022 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00023 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00024 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00025 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00026 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00027 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00028 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00029 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00030 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00031 // 
00032 // This file is part of the PLearn library. For more information on the PLearn
00033 // library, go to the PLearn Web site at www.plearn.org
00034 
00035 /* *******************************************************      
00036  * $Id: BallTreeNearestNeighbors.cc 6861 2007-04-09 19:04:15Z saintmlx $ 
00037  ******************************************************* */
00038 
00039 // Authors: Pascal Lamblin & Marius Muja
00040 
00043 #include "BallTreeNearestNeighbors.h"
00044 #include <plearn/base/lexical_cast.h>
00045 
00046 namespace PLearn {
00047 using namespace std;
00048 
00049 BallTreeNearestNeighbors::BallTreeNearestNeighbors() 
00050     : rmin( 1 ),
00051       train_method( "anchor" )
00052 {
00053     num_neighbors = 1;
00054     expdir = "";
00055     stage = 0;
00056     nstages = -1;
00057     report_progress = 0;
00058 }
00059 
00060 BallTreeNearestNeighbors::BallTreeNearestNeighbors( const VMat& tr_set, const BinBallTree& b_tree )
00061     : rmin( 1 ),
00062       train_method( "anchor" )
00063 {
00064     num_neighbors = 1;
00065     expdir = "";
00066     stage = 1;
00067     nstages = 1;
00068     report_progress = 0;
00069 
00070     setTrainingSet( tr_set );
00071     ball_tree = b_tree;
00072 }
00073 
00074 PLEARN_IMPLEMENT_OBJECT( BallTreeNearestNeighbors, 
00075                          "Organizes hierarchically a set of points to perform efficient  KNN search", 
00076                          "This learner builds a Ball Tree, a hierarchized structure\n"
00077                          "allowing to perform efficient KNN search.\n"
00078                          "Output is formatted as in GenericNearestNeighbors.\n"
00079                          "The square distance to this point can be computed as the error.\n" );
00080 
00081 void BallTreeNearestNeighbors::declareOptions( OptionList& ol )
00082 {
00083     // build options
00084     declareOption( ol, "point_indices", &BallTreeNearestNeighbors::point_indices, 
00085                    OptionBase::buildoption,
00086                    "Indices of the points we will consider" );
00087 
00088     declareOption( ol, "rmin", &BallTreeNearestNeighbors::rmin, OptionBase::buildoption,
00089                    "Max number of points in a leaf node of the tree" );
00090 
00091     declareOption( ol, "train_method", &BallTreeNearestNeighbors::train_method, 
00092                    OptionBase::buildoption,
00093                    "Method used to build the tree. Just one is supported:\n"
00094                    "  \"anchor\" (middle-out building based on Anchor\'s hierarchy\n"
00095         );
00096 
00097     declareOption( ol, "anchor_set", &BallTreeNearestNeighbors::anchor_set, 
00098                    OptionBase::learntoption, 
00099                    "Set of anchors, hierarchizing the set of points" );
00100 
00101     declareOption( ol, "pivot_indices", &BallTreeNearestNeighbors::pivot_indices, 
00102                    OptionBase::learntoption, "Indices of the anchors' centers" );
00103 
00104     // saved options
00105     declareOption( ol, "train_set", &BallTreeNearestNeighbors::train_set, 
00106                    OptionBase::buildoption,
00107                    "Indexed set of points we will be working with" );
00108 
00109     declareOption( ol, "nb_train_points", &BallTreeNearestNeighbors::nb_train_points, 
00110                    OptionBase::learntoption, "Number of points in train_set" );
00111 
00112     declareOption( ol, "nb_points", &BallTreeNearestNeighbors::nb_points, 
00113                    OptionBase::learntoption, "Number of points in point_indices" );
00114 
00115     declareOption( ol, "ball_tree", &BallTreeNearestNeighbors::ball_tree, 
00116                    OptionBase::learntoption, "Built ball-tree" );
00117 
00118 
00119     // Now call the parent class' declareOptions
00120     inherited::declareOptions( ol );
00121 }
00122 
00123 void BallTreeNearestNeighbors::build_()
00124 {
00125     if (train_set) {
00126         // initialize nb_train_points
00127         nb_train_points = train_set.length();
00128         
00129         // if point_indices isn't specified, we take all the points in train_set
00130         if( !point_indices )
00131             point_indices = TVec<int>( 0, nb_train_points-1, 1 );
00132 
00133         // initialize nb_points
00134         nb_points = point_indices.size();
00135     }
00136 }
00137 
00138 
00139 void BallTreeNearestNeighbors::build()
00140 {
00141     inherited::build();
00142     build_();
00143 }
00144 
00145 
00146 void BallTreeNearestNeighbors::makeDeepCopyFromShallowCopy(CopiesMap& copies)
00147 {
00148     inherited::makeDeepCopyFromShallowCopy(copies);
00149 
00150     deepCopyField( ball_tree, copies );
00151     deepCopyField( point_indices, copies );
00152     deepCopyField( anchor_set, copies );
00153     deepCopyField( pivot_indices, copies );
00154 }
00155 
00156 
00157 void BallTreeNearestNeighbors::forget()
00158 {
00161 
00162     anchor_set.resize( 0 );
00163     pivot_indices.resize( 0 );
00164     ball_tree = new BinaryBallTree;
00165     stage = 0;
00166     build();
00167 }
00168 
00169 
00170 
00171 
00172 void BallTreeNearestNeighbors::train()
00173 {
00174     // The role of the train method is to bring the learner up to stage==nstages,
00175     // updating train_stats with training costs measured on-line in the process.
00176 
00177     if( train_method == "anchor" )
00178     {
00179         anchorTrain();
00180     }
00181     else
00182         PLERROR( "train_method \"%s\" not implemented", train_method.c_str() );
00183 }
00184 
00185 
00186 void BallTreeNearestNeighbors::anchorTrain()
00187 {
00188     /*  nstages and stage conventions, for "anchor" train method:
00189      *
00190      *  nstages == -1
00191      *    We will construct ball_tree recursively,
00192      *    until, for all leaf, nb_points <= rmin,
00193      *    no matter how many iterations it will take.
00194      *
00195      *  nstages == 0
00196      *    We want the PLearner il its fresh, blank state.
00197      *
00198      *  nstages == 1
00199      *    We want ball_tree to be a unique leaf node,
00200      *    containing all the point indices, with no children.
00201      *
00202      *  nstages > 1
00203      *    We want to build ball_tree recursively,
00204      *    but limiting the levels of recursion.
00205      *    This means we will decrement this number at each recursive call,
00206      *    the recursion will stop when nstages == 1 or nb_points <= rmin.
00207      *
00208      *  stage == 0
00209      *    The learner is it its fresh, blank state.
00210      *
00211      *  stage == 1
00212      *    The learner has one anchor, that's all.
00213      *
00214      *  Other values of stage might be used one day or anoter...
00215      */
00216 
00217     if( stage == 0 && nstages !=0 )
00218     {
00219         // That means we weren't provided with any anchor nor node parameter,
00220         // or that they were just bullsh!t
00221 
00222         // So, we build a single anchor
00223         pivot_indices.resize( 1 );
00224         pivot_indices[ 0 ] = 0;
00225         Vec pivot = train_set.getSubRow( 0, inputsize() );
00226 
00227         distance_kernel->setDataForKernelMatrix( train_set );
00228         distance_kernel->build();
00229         Vec distances_from_pivot( nb_train_points );
00230         distance_kernel->evaluate_all_i_x( pivot, distances_from_pivot );
00231 
00232         anchor_set.resize( 1 );
00233         Mat* p_anchor = &anchor_set[ 0 ];
00234         p_anchor->resize( nb_points, 2 );
00235         p_anchor->column( 0 ) << Vec( 0, nb_points-1, 1 );
00236         p_anchor->column( 1 ) << distances_from_pivot;
00237         sortRows( *p_anchor, TVec<int>( 1, 1 ), false );
00238 
00239         // then, we build the corresponding tree
00240         ball_tree = leafFromAnchor( 0 );
00241 
00242         ++stage;
00243     }
00244 
00245     if( nstages == 0 )
00246     {
00247         // we want a fresh, blank learner
00248         forget();
00249     }
00250     else if( nstages == 1 )
00251     {
00252         // We have an anchor, and we want a leaf node
00253         ball_tree = leafFromAnchor( 0 );
00254     }
00255     else
00256     {
00257         // nstages to be used on children learners
00258         int new_nstages = nstages<0 ? -1 : nstages-1;
00259 
00260         // First create sqrt( R )-1 anchors, from the initial anchor_set
00261         int nb_anchors = (int) sqrt( (float) nb_points ) + 1 ;
00262         nb_anchors = min( nb_anchors, nb_points );
00263 
00264         createAnchors( nb_anchors-1 ); // because we already have one
00265 
00266         // Convert them into leaf nodes
00267         TVec< BinBallTree > leaf_set = TVec<BinBallTree>( nb_anchors );
00268         for ( int i=0 ; i<nb_anchors ; i++ )
00269         {
00270             leaf_set[ i ] = leafFromAnchor( i );
00271         }
00272 
00273         // Then, group them to form the ball_tree
00274         // keep an index of the leaves
00275         ball_tree = treeFromLeaves( leaf_set );
00276 
00277         // Now, recurse...
00278         for( int i=0 ; i<leaf_set.size() ; i++ )
00279         {
00280             int rec_nb_points = anchor_set[ i ].length();
00281 
00282             // if the leaf is too small, don't do anything
00283             if( rec_nb_points > rmin )
00284             {
00285                 // child learner
00286                 PP<BallTreeNearestNeighbors> p_rec_learner = new BallTreeNearestNeighbors();
00287 
00288                 // initializes child's nstages (see explanation above)
00289                 stringstream out;
00290                 out << new_nstages;
00291                 p_rec_learner->setOption( "nstages" , out.str() );
00292 
00293                 // keep the same training set: it give us all the point coordinates !
00294                 // but we don't want to call forget() after that
00295                 p_rec_learner->setTrainingSet( train_set, false );
00296 
00297                 // however, we only work on the points contained by current leaf
00298                 p_rec_learner->anchor_set.resize( 1 );
00299                 p_rec_learner->anchor_set[ 0 ].resize( rec_nb_points, 2 );
00300                 p_rec_learner->anchor_set[ 0 ] << anchor_set[ i ];
00301 
00302                 p_rec_learner->pivot_indices.resize( 1 );
00303                 p_rec_learner->pivot_indices[ 0 ] = pivot_indices[ i ];
00304 
00305                 p_rec_learner->point_indices.resize( rec_nb_points );
00306                 p_rec_learner->point_indices << 
00307                     p_rec_learner->anchor_set[ 0 ].column( 0 );
00308 
00309                 p_rec_learner->stage = 1; 
00310                 // faudra peut-etre faire ça plus subtilement
00311 
00312                 p_rec_learner->rmin = rmin;
00313                 p_rec_learner->train_method = train_method;
00314                 p_rec_learner->build();
00315                 p_rec_learner->train();
00316 
00317                 // once the child learner is trained, we can get the sub-tree,
00318                 // and link it correctly
00319                 BinBallTree subtree = p_rec_learner->getBallTree();
00320                 leaf_set[ i ]->pivot = subtree->pivot;
00321                 leaf_set[ i ]->radius = subtree->radius;
00322                 leaf_set[ i ]->point_set.resize( subtree->point_set.size() );
00323                 leaf_set[ i ]->point_set << subtree->point_set;
00324                 leaf_set[ i ]->setFirstChild( subtree->getFirstChild() );
00325                 leaf_set[ i ]->setSecondChild( subtree->getSecondChild() );
00326 
00327             }
00328         }
00329     }
00330 }
00331 
00332 
00333 void BallTreeNearestNeighbors::createAnchors( int nb_anchors )
00334 {
00335     // This method creates nb_anchors new anchors, and adds them to anchor_set
00336 
00337     // Make room
00338     int anchor_set_size = anchor_set.size();
00339     anchor_set.resize( anchor_set_size, nb_anchors );
00340 
00341     for( int i=0 ; i<nb_anchors ; i++ )
00342     {
00343         Mat new_anchor = Mat( 1, 2 );
00344         int new_pivot_index;
00345 
00346         // Search for the largest ball.
00347         // pivot of the new anchor will be the point of this ball
00348         // that is the furthest from the pivot.
00349         int largest_index = 0;
00350         real largest_radius = 0;
00351         for( int j=0 ; j<anchor_set_size ; j++ )
00352         {
00353             // points are sorted in decreasing order of distance, 
00354             // so anchor_set[ j ]( 0, 1 ) is the furthest point from 
00355             // pivot_indices[ j ]
00356             real current_radius = anchor_set[ j ]( 0, 1 );
00357             if( current_radius > largest_radius )
00358             {
00359                 largest_radius = current_radius;
00360                 largest_index = j;
00361             }
00362         }
00363 
00364         Mat* p_largest_anchor = &anchor_set[ largest_index ];
00365         new_pivot_index = (int) (*p_largest_anchor)( 0, 0 );
00366 
00367         // assign the point to its new anchor
00368         new_anchor( 0, 0 ) = new_pivot_index;
00369         new_anchor( 0, 1 ) = 0;
00370         Vec new_pivot = train_set.getSubRow( new_pivot_index, inputsize() );
00371 
00372         int largest_anchor_length = p_largest_anchor->length();
00373 
00374         // Verify that largest_anchor owns at least 2 points
00375         if( largest_anchor_length <= 1 )
00376         {
00377             PLERROR("In BallTreeNearestNeighbors::createAnchors, more anchors asked than points");
00378         }
00379 
00380         // delete this point from its original anchor
00381         *p_largest_anchor = p_largest_anchor->
00382             subMatRows( 1, largest_anchor_length-1 );
00383 
00384         // now, try to steal points from all the existing anchors
00385         for( int j=0 ; j<anchor_set_size ; j++ )
00386         {
00387             Mat* p_anchor = &anchor_set[ j ];
00388             int nb_points = p_anchor->length();
00389             int pivot_index = pivot_indices[ j ];
00390             Vec pivot = train_set.getSubRow( pivot_index, inputsize() );
00391             real pivot_pow_dist = powdistance( new_pivot, pivot, 2 );
00392 
00393             // loop on the anchor's points
00394             for( int k=0 ; k<nb_points ; k++ )
00395             {
00396                 int point_index = (int) (*p_anchor)( k, 0 );
00397                 real point_pow_dist = (*p_anchor)( k, 1 );
00398 
00399                 // if this inequality is verified,
00400                 // then we're sure that all the points closer to the pivot 
00401                 // belong to the pivot, and we don't need to check
00402                 if( 4*point_pow_dist < pivot_pow_dist )
00403                 {
00404                     break;
00405                 }
00406 
00407                 Vec point = train_set.getSubRow( point_index, inputsize() );
00408                 real new_pow_dist = powdistance( new_pivot, point, 2 );
00409 
00410                 // if the point is closer to the new pivot, then steal it
00411                 if( new_pow_dist < point_pow_dist )
00412                 {
00413                     Vec new_row( 2 );
00414                     new_row[ 0 ] = point_index;
00415                     new_row[ 1 ] = new_pow_dist;
00416                     new_anchor.appendRow( new_row );
00417 
00418                     *p_anchor = removeRow( *p_anchor, k );
00419                     // bleaah, this is ugly !
00420                     --k;
00421                     --nb_points;
00422                 }
00423             }
00424         }
00425 
00426         // sort the points by decreasing distance
00427         sortRows( new_anchor, TVec<int>( 1, 1 ), false );
00428 
00429         // append the new anchor to the anchor_set (and same for pivot)
00430         anchor_set.append( new_anchor );
00431         pivot_indices.append( new_pivot_index );
00432         ++anchor_set_size;
00433     }
00434 }
00435 
00436 BinBallTree BallTreeNearestNeighbors::leafFromAnchor( int anchor_index )
00437 {
00438     BinBallTree leaf = new BinaryBallTree();
00439 
00440     int pivot_index = pivot_indices[ anchor_index ];
00441     leaf->pivot = train_set.getSubRow( pivot_index, inputsize() );
00442 
00443     leaf->radius = anchor_set[ anchor_index ]( 0, 1 );
00444 
00445     int nb_leaf_points = anchor_set[ anchor_index ].length();
00446     leaf->point_set.resize( nb_leaf_points );
00447     leaf->point_set << anchor_set[ anchor_index ].column( 0 );
00448 
00449     return leaf;
00450 }
00451 
00452 
00453 BinBallTree BallTreeNearestNeighbors::treeFromLeaves( const TVec<BinBallTree>& leaves )
00454 {
00455     int nb_nodes = leaves.size();
00456     TVec<BinBallTree> nodes = TVec<BinBallTree>( nb_nodes );
00457     nodes << leaves;
00458 
00459     // if there is no leaf
00460     if( nb_nodes < 1 )
00461     {
00462         PLERROR( "In BallTreeNearestNeighbors::treeFromLeaves(): no leaf existing" );
00463     }
00464 
00465     while( nb_nodes > 1 )
00466     {
00467         int min_i = 0;
00468         int min_j = 0;
00469         Vec min_center;
00470         real min_radius = -1;
00471 
00472         // we get the most "compatible" pair of nodes :
00473         // the ball containing them both is the smallest
00474         for( int i=0 ; i<nb_nodes ; i++ )
00475         {
00476             Vec center_i = nodes[ i ]->pivot;
00477             real radius_i = nodes[ i ]->radius;
00478 
00479             // to scan all pairs only once, and avoid i==j
00480             for( int j=0 ; j<i ; j++ )
00481             {
00482                 Vec center_j = nodes[ j ]->pivot;
00483                 real radius_j = nodes[ j ]->radius;
00484 
00485                 Vec t_center;
00486                 real t_radius;
00487                 smallestContainer( center_i, radius_i, center_j, radius_j, 
00488                                    t_center, t_radius );
00489 
00490                 if( t_radius < min_radius || min_radius < 0 )
00491                 {
00492                     min_i = i;
00493                     min_j = j ;
00494                     min_radius = t_radius;
00495                     min_center = t_center;
00496                 }
00497             }
00498         }
00499 
00500 #ifdef DEBUG_CHECK_NAN
00501         if (min_center.hasMissing())
00502             PLERROR("In BallTreeNearestNeighbors::treeFromLeaves: min_center is NaN");
00503 #endif
00504         
00505         // Group these two nodes into a parent_node.
00506         // TODO: something more sensible for the radius and center...
00507         BinBallTree parent_node = new BinaryBallTree();
00508         parent_node->pivot = min_center;
00509         parent_node->radius = min_radius;
00510         parent_node->setFirstChild( nodes[ min_i ] );
00511         parent_node->setSecondChild( nodes[ min_j ] );
00512 
00513         nodes[ min_j ] = parent_node;
00514         nodes.remove( min_i );
00515 
00516         --nb_nodes;
00517     }
00518 
00519     // then, we have only one anchor
00520     BinBallTree root = nodes[ 0 ];
00521     return root;
00522 }
00523 
00524 
00525 BinBallTree BallTreeNearestNeighbors::getBallTree()
00526 {
00527     return ball_tree;
00528 }
00529 
00530 
00531 void BallTreeNearestNeighbors::computeOutputAndCosts(
00532     const Vec& input, const Vec& target, Vec& output, Vec& costs ) const
00533 {
00534     int nout = outputsize();
00535     output.resize( nout );
00536     costs.resize( num_neighbors );
00537 
00538     // we launch a k-nearest-neighbors query on the root node (ball_tree)
00539     priority_queue< pair<real,int> > q;
00540     FindBallKNN( q, input, num_neighbors );
00541 
00542     // dequeue the found nearest neighbors, beginning by the farthest away
00543     int n_found = int(q.size());
00544     TVec<int> neighbors( n_found );
00545     for( int i=n_found-1 ; i>=0 ; i-- )
00546     {
00547         const pair<real,int>& cur_top = q.top();
00548         costs[i] = cur_top.first;
00549         neighbors[i] = cur_top.second;
00550         q.pop();
00551     }
00552 
00553     // fill costs with missing values
00554     for( int i= n_found ; i<num_neighbors ; i++ )
00555         costs[i] = MISSING_VALUE;
00556 
00557     constructOutputVector( neighbors, output );
00558 }
00559 
00560 void BallTreeNearestNeighbors::computeOutput(
00561     const Vec& input, Vec& output ) const
00562 {
00563     // Compute the output from the input.
00564     // int nout = outputsize();
00565     // output.resize(nout);
00566 
00567     int nout = outputsize();
00568     output.resize( nout );
00569 
00570     // we launch a k-nearest-neighbors query on the root node (ball_tree)
00571     priority_queue< pair<real,int> > q;
00572     FindBallKNN( q, input, num_neighbors );
00573 
00574     // dequeue the found nearest neighbors, beginning by the farthest away
00575     int n_found = int(q.size());
00576     TVec<int> neighbors( n_found );
00577     for( int i=n_found-1 ; i>=0 ; i-- )
00578     {
00579         const pair<real,int>& cur_top = q.top();
00580         neighbors[i] = cur_top.second;
00581         q.pop();
00582     }
00583 
00584     constructOutputVector( neighbors, output );
00585 
00586 }
00587 
00588 
00589 void BallTreeNearestNeighbors::computeCostsFromOutputs(
00590     const Vec& input, const Vec& output, const Vec& target, Vec& costs ) const
00591 {
00592     // Compute the costs from *already* computed output.
00593     costs.resize( num_neighbors );
00594 
00595     int inputsize = train_set->inputsize();
00596     int targetsize = train_set->targetsize();
00597     int weightsize = train_set->weightsize();
00598 
00599     Mat out( num_neighbors, inputsize );
00600 
00601     if( copy_input )
00602     {
00603         for( int i=0 ; i<num_neighbors ; i++ )
00604             out( i ) << output.subVec( i*outputsize(), inputsize );
00605     }
00606     else if( copy_index )
00607     {
00608         int offset = 0;
00609 
00610         if( copy_target )
00611             offset += targetsize;
00612 
00613         if( copy_weight )
00614             offset += weightsize;
00615 
00616         for( int i=0 ; i<num_neighbors ; i++ )
00617             out( i ) << train_set( (int) output[ i*outputsize() + offset ] );
00618     }
00619     else
00620     {
00621         PLERROR( "computeCostsFromOutput:\n"
00622                  "neither indices nor coordinates of output computed\n" );
00623     }
00624 
00625     for( int i=0 ; i<num_neighbors ; i++ )
00626         costs[ i ] = powdistance( input, out( i ) );
00627 }
00628 
00629 TVec<string> BallTreeNearestNeighbors::getTestCostNames() const
00630 {
00631     return TVec<string>( num_neighbors, "squared_distance" );
00632 }
00633 
00634 TVec<string> BallTreeNearestNeighbors::getTrainCostNames() const
00635 {
00636     return TVec<string>();
00637 }
00638 
00639 bool BallTreeNearestNeighbors::intersect(
00640     const Vec& center1, const real& powrad1,
00641     const Vec& center2, const real& powrad2 )
00642 {
00643     real radius1 = sqrt( powrad1 );
00644     real radius2 = sqrt( powrad2 );
00645 
00646     real pow_dist = powdistance( center1, center2, 2 );
00647     real rad_sum = radius1 + radius2;
00648     bool result = ( pow_dist <= ( rad_sum * rad_sum ) );
00649     return result;
00650 }
00651 
00652 bool BallTreeNearestNeighbors::contain(
00653     const Vec& center1, const real& powrad1,
00654     const Vec& center2, const real& powrad2 )
00655 {
00656     real radius1 = sqrt( powrad1 );
00657     real radius2 = sqrt( powrad2 );
00658     real rad_dif = radius1 - radius2;
00659 
00660     if( rad_dif >= 0 )
00661     {
00662         real pow_dist = powdistance( center1, center2, 2 );
00663         bool result = ( pow_dist <= ( rad_dif * rad_dif ) );
00664         return result;
00665     }
00666     else
00667     {
00668         return false;
00669     }
00670 }
00671 
00672 void BallTreeNearestNeighbors::smallestContainer(
00673     const Vec& center1, const real& powrad1,
00674     const Vec& center2, const real& powrad2,
00675     Vec& t_center, real& t_powrad )
00676 {
00677     if( center1 == center2 )
00678     {
00679         t_center = center1;
00680         t_powrad = max( powrad1, powrad2 );
00681     }
00682     else if( contain( center1, powrad1, center2, powrad2 ) )
00683     {
00684         t_center = center1;
00685         t_powrad = powrad1;
00686     }
00687     else if( contain( center2, powrad2, center1, powrad1 ) )
00688     {
00689         t_center = center2;
00690         t_powrad = powrad2;
00691     }
00692     else
00693     {
00694         real radius1 = sqrt( powrad1 );
00695         real radius2 = sqrt( powrad2 );
00696         real center_dist = dist( center1, center2, 2 ) ;
00697         real coef = ( radius1 - radius2 ) / center_dist ;
00698         t_center = real(0.5) * ( ( 1 + coef ) * center1  +  ( 1 - coef ) * center2 ) ;
00699         real t_radius = real(0.5) * ( center_dist + radius1 + radius2 ) ;
00700         t_powrad = t_radius * t_radius;
00701     }
00702 
00703 #ifdef DEBUG_CHECK_NAN
00704     if (t_center.hasMissing())
00705         PLERROR("In BallTreeNearestNeighbors::smallestContainer: t_center is NaN.");
00706 #endif
00707 }
00708 
00709 
00710 
00711 void BallTreeNearestNeighbors::BallKNN(
00712      priority_queue< pair<real,int> >& q, BinBallTree node,
00713      const Vec& t, real& d2_sofar, real d2_pivot, const int k ) const
00714 {
00715     real d_minp = max( sqrt(d2_pivot) - node->radius, real(0.0) );
00716 #ifdef DEBUG_CHECK_NAN
00717     if (isnan(d_minp))
00718         PLERROR("BallTreeNearestNeighbors::BallKNN: d_minp is NaN");
00719 #endif
00720 
00721     if (d_minp*d_minp > d2_sofar)
00722     {
00723         // no chance of finding anything closer around this node
00724         return;
00725     }
00726     else if (node->point_set.size()!=0) // node is leaf
00727     {
00728         int n_points = node->point_set.size();
00729         for( int i=0 ; i<n_points ; i++ )
00730         {
00731             int j = node->point_set[i];
00732             real dist;
00733             // last point is pivot, and we already now the distance
00734             if( i==n_points-1 )
00735             {
00736                 dist = d2_pivot;
00737             }
00738             else
00739             {
00740                 Vec x = train_set.getSubRow(j, inputsize());
00741                 dist = powdistance(x, t, 2);
00742             }
00743             if( dist < d2_sofar )
00744             {
00745                 q.push( make_pair(dist, j) );
00746                 int n_found = int(q.size());
00747                 if( n_found > k )
00748                     q.pop();
00749                 if( n_found >= k )
00750                     d2_sofar = q.top().first;
00751             }
00752         }
00753     }
00754     else if (!node->isEmpty()) // node is not leaf
00755     {
00756         BinBallTree node1 = node->getFirstChild();
00757         BinBallTree node2 = node->getSecondChild();
00758 
00759         real d2_pivot1 = powdistance(t, node1->pivot, 2);
00760         real d2_pivot2 = powdistance(t, node2->pivot, 2);
00761 
00762         if( d2_pivot1 > d2_pivot2 ) // node1 is closer to t
00763         {
00764             pl_swap(node1, node2);
00765             pl_swap(d2_pivot1, d2_pivot2);
00766         }
00767 
00768         BallKNN(q, node1, t, d2_sofar, d2_pivot1, k);
00769         BallKNN(q, node2, t, d2_sofar, d2_pivot2, k); 
00770     }
00771 }
00772 
00773 
00774 void BallTreeNearestNeighbors::FindBallKNN(
00775     priority_queue< pair<real,int> >& q, const Vec& point, const int k ) const
00776 {
00777     real d2_sofar;
00778     pl_isnumber("+inf", &d2_sofar);
00779     real d2_pivot = powdistance(point, ball_tree->pivot, 2);
00780 //    real d_minp = 0;
00781     BallKNN(q, ball_tree, point, d2_sofar, d2_pivot, k);
00782 }
00783 
00784 } // end of namespace PLearn
00785 
00786 
00787 /*
00788   Local Variables:
00789   mode:c++
00790   c-basic-offset:4
00791   c-file-style:"stroustrup"
00792   c-file-offsets:((innamespace . 0)(inline-open . 0))
00793   indent-tabs-mode:nil
00794   fill-column:79
00795   End:
00796 */
00797 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines