FastJet 3.0.2
ClosestPair2D.cc
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends