Main Page | Namespace List | Class Hierarchy | Class List | Directories | File List | Namespace Members | Class Members | File Members

ClosestPair2D.cc

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

Generated on Mon Apr 2 20:57:48 2007 for fastjet by  doxygen 1.4.2