PLearn 0.1
|
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 :