FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: ClosestPair2D.cc 2577 2011-09-13 15:11:38Z salam $ 00003 // 00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 00005 // 00006 //---------------------------------------------------------------------- 00007 // This file is part of FastJet. 00008 // 00009 // FastJet is free software; you can redistribute it and/or modify 00010 // it under the terms of the GNU General Public License as published by 00011 // the Free Software Foundation; either version 2 of the License, or 00012 // (at your option) any later version. 00013 // 00014 // The algorithms that underlie FastJet have required considerable 00015 // development and are described in hep-ph/0512210. If you use 00016 // FastJet as part of work towards a scientific publication, please 00017 // include a citation to the FastJet paper. 00018 // 00019 // FastJet is distributed in the hope that it will be useful, 00020 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00021 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00022 // GNU General Public License for more details. 00023 // 00024 // You should have received a copy of the GNU General Public License 00025 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00026 //---------------------------------------------------------------------- 00027 //ENDHEADER 00028 00029 #include "fastjet/internal/ClosestPair2D.hh" 00030 00031 #include<limits> 00032 #include<iostream> 00033 #include<iomanip> 00034 #include<algorithm> 00035 00036 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00037 00038 const unsigned int huge_unsigned = 4294967295U; 00039 const unsigned int twopow31 = 2147483648U; 00040 00041 using namespace std; 00042 00043 //---------------------------------------------------------------------- 00044 /// takes a point and sets a shuffle with the given shift... 00045 void ClosestPair2D::_point2shuffle(Point & point, Shuffle & shuffle, 00046 unsigned int shift) { 00047 00048 Coord2D renorm_point = (point.coord - _left_corner)/_range; 00049 // make sure the point is sensible 00050 //cerr << point.coord.x <<" "<<point.coord.y<<endl; 00051 assert(renorm_point.x >=0); 00052 assert(renorm_point.x <=1); 00053 assert(renorm_point.y >=0); 00054 assert(renorm_point.y <=1); 00055 00056 shuffle.x = static_cast<unsigned int>(twopow31 * renorm_point.x) + shift; 00057 shuffle.y = static_cast<unsigned int>(twopow31 * renorm_point.y) + shift; 00058 shuffle.point = &point; 00059 } 00060 00061 //---------------------------------------------------------------------- 00062 /// compares this shuffle with the other one 00063 bool ClosestPair2D::Shuffle::operator<(const Shuffle & q) const { 00064 00065 if (floor_ln2_less(x ^ q.x, y ^ q.y)) { 00066 // i = 2 in Chan's algorithm 00067 return (y < q.y); 00068 } else { 00069 // i = 1 in Chan's algorithm 00070 return (x < q.x); 00071 } 00072 } 00073 00074 00075 00076 //---------------------------------------------------------------------- 00077 void ClosestPair2D::_initialize(const std::vector<Coord2D> & positions, 00078 const Coord2D & left_corner, 00079 const Coord2D & right_corner, 00080 unsigned int max_size) { 00081 00082 unsigned int n_positions = positions.size(); 00083 assert(max_size >= n_positions); 00084 00085 //_points(positions.size()) 00086 00087 // allow the points array to grow to the following size 00088 _points.resize(max_size); 00089 // currently unused points are immediately made available on the 00090 // stack 00091 for (unsigned int i = n_positions; i < max_size; i++) { 00092 _available_points.push(&(_points[i])); 00093 } 00094 00095 _left_corner = left_corner; 00096 _range = max((right_corner.x - left_corner.x), 00097 (right_corner.y - left_corner.y)); 00098 00099 // initialise the coordinates for the points and create the zero-shifted 00100 // shuffle array 00101 vector<Shuffle> shuffles(n_positions); 00102 for (unsigned int i = 0; i < n_positions; i++) { 00103 // set up the points 00104 _points[i].coord = positions[i]; 00105 _points[i].neighbour_dist2 = numeric_limits<double>::max(); 00106 _points[i].review_flag = 0; 00107 00108 // create shuffle with 0 shift. 00109 _point2shuffle(_points[i], shuffles[i], 0); 00110 } 00111 00112 // establish what our shifts will be 00113 for (unsigned ishift = 0; ishift < _nshift; ishift++) { 00114 // make sure we use double-precision for calculating the shifts 00115 // since otherwise we will hit integer overflow. 00116 _shifts[ishift] = static_cast<unsigned int>(((twopow31*1.0)*ishift)/_nshift); 00117 if (ishift == 0) {_rel_shifts[ishift] = 0;} 00118 else {_rel_shifts[ishift] = _shifts[ishift] - _shifts[ishift-1];} 00119 } 00120 //_shifts[0] = 0; 00121 //_shifts[1] = static_cast<unsigned int>((twopow31*1.0)/3.0); 00122 //_shifts[2] = static_cast<unsigned int>((twopow31*2.0)/3.0); 00123 //_rel_shifts[0] = 0; 00124 //_rel_shifts[1] = _shifts[1]; 00125 //_rel_shifts[2] = _shifts[2]-_shifts[1]; 00126 00127 // and how we will search... 00128 //_cp_search_range = 49; 00129 _cp_search_range = 30; 00130 _points_under_review.reserve(_nshift * _cp_search_range); 00131 00132 // now initialise the three trees 00133 for (unsigned int ishift = 0; ishift < _nshift; ishift++) { 00134 00135 // shift the shuffles if need be. 00136 if (ishift > 0) { 00137 unsigned rel_shift = _rel_shifts[ishift]; 00138 for (unsigned int i = 0; i < shuffles.size(); i++) { 00139 shuffles[i] += rel_shift; } 00140 } 00141 00142 // sort the shuffles 00143 sort(shuffles.begin(), shuffles.end()); 00144 00145 // and create the search tree 00146 _trees[ishift] = auto_ptr<Tree>(new Tree(shuffles, max_size)); 00147 00148 // now we look for the closest-pair candidates on this tree 00149 circulator circ = _trees[ishift]->somewhere(), start=circ; 00150 // the actual range in which we search 00151 unsigned int CP_range = min(_cp_search_range, n_positions-1); 00152 do { 00153 Point * this_point = circ->point; 00154 //cout << _ID(this_point) << " "; 00155 this_point->circ[ishift] = circ; 00156 // run over all points within _cp_search_range of this_point on tree 00157 circulator other = circ; 00158 for (unsigned i=0; i < CP_range; i++) { 00159 ++other; 00160 double dist2 = this_point->distance2(*other->point); 00161 if (dist2 < this_point->neighbour_dist2) { 00162 this_point->neighbour_dist2 = dist2; 00163 this_point->neighbour = other->point; 00164 } 00165 } 00166 } while (++circ != start); 00167 //cout << endl<<endl; 00168 } 00169 00170 // now initialise the heap object... 00171 vector<double> mindists2(n_positions); 00172 for (unsigned int i = 0; i < n_positions; i++) { 00173 mindists2[i] = _points[i].neighbour_dist2;} 00174 00175 _heap = auto_ptr<MinHeap>(new MinHeap(mindists2, max_size)); 00176 } 00177 00178 00179 //----------------------------------------------------------------------= 00180 void ClosestPair2D::closest_pair(unsigned int & ID1, unsigned int & ID2, 00181 double & distance2) const { 00182 ID1 = _heap->minloc(); 00183 ID2 = _ID(_points[ID1].neighbour); 00184 distance2 = _points[ID1].neighbour_dist2; 00185 if (ID1 > ID2) swap(ID1,ID2); 00186 } 00187 00188 00189 //---------------------------------------------------------------------- 00190 inline void ClosestPair2D::_add_label(Point * point, unsigned int review_flag) { 00191 00192 // if it's not already under review, then put it on the list of 00193 // points needing review 00194 if (point->review_flag == 0) _points_under_review.push_back(point); 00195 00196 // OR the point's current flag with the requested review flag 00197 point->review_flag |= review_flag; 00198 } 00199 00200 //---------------------------------------------------------------------- 00201 inline void ClosestPair2D::_set_label(Point * point, unsigned int review_flag) { 00202 00203 // if it's not already under review, then put it on the list of 00204 // points needing review 00205 if (point->review_flag == 0) _points_under_review.push_back(point); 00206 00207 // SET the point's current flag to the requested review flag 00208 point->review_flag = review_flag; 00209 } 00210 00211 //---------------------------------------------------------------------- 00212 void ClosestPair2D::remove(unsigned int ID) { 00213 00214 //cout << "While removing " << ID <<endl; 00215 00216 Point * point_to_remove = & (_points[ID]); 00217 00218 // remove this point from the search tree 00219 _remove_from_search_tree(point_to_remove); 00220 00221 // the above statement labels certain points as needing "review" -- 00222 // deal with them... 00223 _deal_with_points_to_review(); 00224 } 00225 00226 00227 //---------------------------------------------------------------------- 00228 void ClosestPair2D::_remove_from_search_tree(Point * point_to_remove) { 00229 00230 // add this point to the list of "available" points (this is 00231 // relevant also from the point of view of determining the range 00232 // over which we circulate). 00233 _available_points.push(point_to_remove); 00234 00235 // label it so that it goes from the heap 00236 _set_label(point_to_remove, _remove_heap_entry); 00237 00238 // establish the range over which we search (a) for points that have 00239 // acquired a new neighbour and (b) for points which had ID as their 00240 // neighbour; 00241 00242 unsigned int CP_range = min(_cp_search_range, size()-1); 00243 00244 // then, for each shift 00245 for (unsigned int ishift = 0; ishift < _nshift; ishift++) { 00246 //cout << " ishift = " << ishift <<endl; 00247 // get the circulator for the point we'll remove and its successor 00248 circulator removed_circ = point_to_remove->circ[ishift]; 00249 circulator right_end = removed_circ.next(); 00250 // remove the point 00251 _trees[ishift]->remove(removed_circ); 00252 00253 // next find the point CP_range points to the left 00254 circulator left_end = right_end, orig_right_end = right_end; 00255 for (unsigned int i = 0; i < CP_range; i++) {left_end--;} 00256 00257 if (size()-1 < _cp_search_range) { 00258 // we have a smaller range now than before -- but when seeing who 00259 // could have had ID as a neighbour, we still need to use the old 00260 // range for seeing how far back we search (but new separation between 00261 // points). [cf CCN28-42] 00262 left_end--; right_end--; 00263 } 00264 00265 // and then for each left-end point: establish if the removed 00266 // point was its neighbour [in which case a new neighbour must be 00267 // found], otherwise see if the right-end point is a closer neighbour 00268 do { 00269 Point * left_point = left_end->point; 00270 00271 //cout << " visited " << setw(3)<<_ID(left_point)<<" (its neighbour was "<< setw(3)<< _ID(left_point->neighbour)<<")"<<endl; 00272 00273 if (left_point->neighbour == point_to_remove) { 00274 // we'll deal with it later... 00275 _add_label(left_point, _review_neighbour); 00276 } else { 00277 // check to see if right point has become its closest neighbour 00278 double dist2 = left_point->distance2(*right_end->point); 00279 if (dist2 < left_point->neighbour_dist2) { 00280 left_point->neighbour = right_end->point; 00281 left_point->neighbour_dist2 = dist2; 00282 // NB: (LESSER) REVIEW NEEDED HERE TOO... 00283 _add_label(left_point, _review_heap_entry); 00284 } 00285 } 00286 ++right_end; 00287 } while (++left_end != orig_right_end); 00288 } // ishift... 00289 00290 } 00291 00292 00293 //---------------------------------------------------------------------- 00294 void ClosestPair2D::_deal_with_points_to_review() { 00295 00296 // the range in which we carry out searches for new neighbours on 00297 // the search tree 00298 unsigned int CP_range = min(_cp_search_range, size()-1); 00299 00300 // now deal with the points that are "under review" in some way 00301 // (have lost their neighbour, or need their heap entry updating) 00302 while(_points_under_review.size() > 0) { 00303 // get the point to be considered 00304 Point * this_point = _points_under_review.back(); 00305 // remove it from the list 00306 _points_under_review.pop_back(); 00307 00308 if (this_point->review_flag & _remove_heap_entry) { 00309 // make sure no other flags are on (it wouldn't be consistent?) 00310 assert(!(this_point->review_flag ^ _remove_heap_entry)); 00311 _heap->remove(_ID(this_point)); 00312 } 00313 // check to see if the _review_neighbour flag is on 00314 else { 00315 if (this_point->review_flag & _review_neighbour) { 00316 this_point->neighbour_dist2 = numeric_limits<double>::max(); 00317 // among all three shifts 00318 for (unsigned int ishift = 0; ishift < _nshift; ishift++) { 00319 circulator other = this_point->circ[ishift]; 00320 // among points within CP_range 00321 for (unsigned i=0; i < CP_range; i++) { 00322 ++other; 00323 double dist2 = this_point->distance2(*other->point); 00324 if (dist2 < this_point->neighbour_dist2) { 00325 this_point->neighbour_dist2 = dist2; 00326 this_point->neighbour = other->point; 00327 } 00328 } 00329 } 00330 } 00331 00332 // for any non-zero review flag we'll have to update the heap 00333 _heap->update(_ID(this_point), this_point->neighbour_dist2); 00334 } 00335 00336 // "delabel" the point 00337 this_point->review_flag = 0; 00338 00339 } 00340 00341 } 00342 00343 //---------------------------------------------------------------------- 00344 unsigned int ClosestPair2D::insert(const Coord2D & new_coord) { 00345 00346 // get hold of a point 00347 assert(_available_points.size() > 0); 00348 Point * new_point = _available_points.top(); 00349 _available_points.pop(); 00350 00351 // set the point's coordinate 00352 new_point->coord = new_coord; 00353 00354 // now find it's neighbour in the search tree 00355 _insert_into_search_tree(new_point); 00356 00357 // sort out other points that may have been affected by this, 00358 // and/or for which the heap needs to be updated 00359 _deal_with_points_to_review(); 00360 00361 // 00362 return _ID(new_point); 00363 } 00364 00365 //---------------------------------------------------------------------- 00366 unsigned int ClosestPair2D::replace(unsigned int ID1, unsigned int ID2, 00367 const Coord2D & position) { 00368 00369 // deletion from tree... 00370 Point * point_to_remove = & (_points[ID1]); 00371 _remove_from_search_tree(point_to_remove); 00372 00373 point_to_remove = & (_points[ID2]); 00374 _remove_from_search_tree(point_to_remove); 00375 00376 // insertion into tree 00377 // get hold of a point 00378 Point * new_point = _available_points.top(); 00379 _available_points.pop(); 00380 00381 // set the point's coordinate 00382 new_point->coord = position; 00383 00384 // now find it's neighbour in the search tree 00385 _insert_into_search_tree(new_point); 00386 00387 // the above statement labels certain points as needing "review" -- 00388 // deal with them... 00389 _deal_with_points_to_review(); 00390 00391 // 00392 return _ID(new_point); 00393 00394 } 00395 00396 00397 //---------------------------------------------------------------------- 00398 void ClosestPair2D::replace_many( 00399 const std::vector<unsigned int> & IDs_to_remove, 00400 const std::vector<Coord2D> & new_positions, 00401 std::vector<unsigned int> & new_IDs) { 00402 00403 // deletion from tree... 00404 for (unsigned int i = 0; i < IDs_to_remove.size(); i++) { 00405 _remove_from_search_tree(& (_points[IDs_to_remove[i]])); 00406 } 00407 00408 // insertion into tree 00409 new_IDs.resize(0); 00410 for (unsigned int i = 0; i < new_positions.size(); i++) { 00411 Point * new_point = _available_points.top(); 00412 _available_points.pop(); 00413 // set the point's coordinate 00414 new_point->coord = new_positions[i]; 00415 // now find it's neighbour in the search tree 00416 _insert_into_search_tree(new_point); 00417 // record the ID 00418 new_IDs.push_back(_ID(new_point)); 00419 } 00420 00421 // the above statement labels certain points as needing "review" -- 00422 // deal with them... 00423 _deal_with_points_to_review(); 00424 00425 } 00426 00427 00428 //---------------------------------------------------------------------- 00429 void ClosestPair2D::_insert_into_search_tree(Point * new_point) { 00430 00431 // this point will have to have it's heap entry reviewed... 00432 _set_label(new_point, _review_heap_entry); 00433 00434 // set the current distance to "infinity" 00435 new_point->neighbour_dist2 = numeric_limits<double>::max(); 00436 00437 // establish how far we will be searching; 00438 unsigned int CP_range = min(_cp_search_range, size()-1); 00439 00440 for (unsigned ishift = 0; ishift < _nshift; ishift++) { 00441 // create the shuffle 00442 Shuffle new_shuffle; 00443 _point2shuffle(*new_point, new_shuffle, _shifts[ishift]); 00444 00445 // insert it into the tree 00446 circulator new_circ = _trees[ishift]->insert(new_shuffle); 00447 new_point->circ[ishift] = new_circ; 00448 00449 // now get hold of the right and left edges of the region we will be 00450 // looking at (cf CCN28-43) 00451 circulator right_edge = new_circ; right_edge++; 00452 circulator left_edge = new_circ; 00453 for (unsigned int i = 0; i < CP_range; i++) {left_edge--;} 00454 00455 // now 00456 do { 00457 Point * left_point = left_edge->point; 00458 Point * right_point = right_edge->point; 00459 00460 // see if the new point is closer to the left-edge than the latter's 00461 // current neighbour 00462 double new_dist2 = left_point->distance2(*new_point); 00463 if (new_dist2 < left_point->neighbour_dist2) { 00464 left_point->neighbour_dist2 = new_dist2; 00465 left_point->neighbour = new_point; 00466 _add_label(left_point, _review_heap_entry); 00467 } 00468 00469 // see if the right-point is closer to the new point than it's current 00470 // neighbour 00471 new_dist2 = new_point->distance2(*right_point); 00472 if (new_dist2 < new_point->neighbour_dist2) { 00473 new_point->neighbour_dist2 = new_dist2; 00474 new_point->neighbour = right_point; 00475 } 00476 00477 // if the right-edge point was the left-edge's neighbour, then 00478 // then it's just gone off-radar and the left-point will need to 00479 // have its neighbour recalculated [actually, this is overdoing 00480 // it a little, since right point may be an less "distant" 00481 // (circulator distance) in one of the other shifts -- but not 00482 // sure how to deal with this...] 00483 if (left_point->neighbour == right_point) { 00484 _add_label(left_point, _review_neighbour); 00485 } 00486 00487 // shift the left and right edges until left edge hits new_circ 00488 right_edge++; 00489 } while (++left_edge != new_circ); 00490 } 00491 } 00492 00493 FASTJET_END_NAMESPACE 00494