FastJet 3.0.2
ClosestPair2D.hh
00001 //STARTHEADER
00002 // $Id: ClosestPair2D.hh 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 #ifndef __FASTJET_CLOSESTPAIR2D__HH__
00030 #define __FASTJET_CLOSESTPAIR2D__HH__
00031 
00032 #include<vector>
00033 #include<stack>
00034 #include<iostream>
00035 #include "fastjet/internal/ClosestPair2DBase.hh"
00036 #include "fastjet/internal/SearchTree.hh"
00037 #include "fastjet/internal/MinHeap.hh"
00038 
00039 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00040 
00041 //----------------------------------------------------------------------
00042 /// \if internal_doc
00043 /// @ingroup internal
00044 /// \class ClosestPair2D
00045 /// concrete implementation for finding closest pairs in 2D -- will
00046 /// use Chan's (hopefully efficient) shuffle based structures
00047 /// \endif
00048 class ClosestPair2D : public ClosestPair2DBase {
00049 public:
00050   /// constructor from a vector of 2D positions -- number of objects
00051   /// after insertion and deletion must never exceed positions.size();
00052   /// objects are given IDs that correspond to their index in the vector 
00053   /// of positions
00054   ClosestPair2D(const std::vector<Coord2D> & positions, 
00055                 const Coord2D & left_corner, const Coord2D & right_corner) {
00056     _initialize(positions, left_corner, right_corner, positions.size());
00057   };
00058 
00059   /// constructor which allows structure to grow beyond positions.size(), up
00060   /// to max_size
00061   ClosestPair2D(const std::vector<Coord2D> & positions, 
00062                 const Coord2D & left_corner, const Coord2D & right_corner,
00063                 const unsigned int max_size) {
00064     _initialize(positions, left_corner, right_corner, max_size);
00065   };
00066 
00067   /// provides the IDs of the closest pair as well as the distance between
00068   /// them
00069   void closest_pair(unsigned int & ID1, unsigned int & ID2, 
00070                     double & distance2) const;
00071 
00072   /// removes the entry labelled by ID from the object;
00073   void remove(unsigned int ID);
00074  
00075   /// inserts the position into the closest pair structure and returns the
00076   /// ID that has been allocated for the object.
00077   unsigned int insert(const Coord2D &);
00078 
00079   /// removes ID1 and ID2 and inserts position, returning the ID 
00080   /// corresponding to position...
00081   virtual unsigned int replace(unsigned int ID1, unsigned int ID2, 
00082                                const Coord2D & position);
00083 
00084   /// replaces IDs_to_remove with points at the new_positions
00085   /// indicating the IDs allocated to the new points in new_IDs
00086   virtual void replace_many(const std::vector<unsigned int> & IDs_to_remove,
00087                             const std::vector<Coord2D> & new_positions,
00088                             std::vector<unsigned int> & new_IDs);
00089 
00090   // mostly for checking how things are working...
00091   inline void print_tree_depths(std::ostream & outdev) const {
00092     outdev    << _trees[0]->max_depth() << " "
00093               << _trees[1]->max_depth() << " "
00094               << _trees[2]->max_depth() << "\n";
00095   };
00096 
00097   unsigned int size();
00098 
00099 private:
00100   
00101   void _initialize(const std::vector<Coord2D> & positions, 
00102               const Coord2D & left_corner, const Coord2D & right_corner,
00103               const unsigned int max_size);
00104 
00105   static const unsigned int _nshift = 3;
00106 
00107   class Point; // will be defined below
00108 
00109   /// since sets of three objects will crop up repeatedly, useful
00110   /// to have a triplet class?
00111   template<class T> class triplet {
00112   public:
00113     inline const T & operator[](unsigned int i) const {return _contents[i];};
00114     inline       T & operator[](unsigned int i)       {return _contents[i];};
00115   private:
00116     T _contents[_nshift];
00117   };
00118 
00119 
00120   /// class that will take care of ordering of shuffles for us
00121   class Shuffle {
00122   public:
00123     unsigned int x, y;
00124     Point * point;
00125     bool operator<(const Shuffle &) const;
00126     void operator+=(unsigned int shift) {x += shift; y+= shift;};
00127   };
00128 
00129   typedef SearchTree<Shuffle>     Tree;
00130   typedef Tree::circulator        circulator;
00131   typedef Tree::const_circulator  const_circulator;
00132 
00133 
00134   triplet<std::auto_ptr<Tree> >  _trees;
00135   std::auto_ptr<MinHeap> _heap;
00136   std::vector<Point>     _points;
00137   std::stack<Point *>    _available_points;
00138 
00139   /// points that are "under review" in some way
00140   std::vector<Point *>   _points_under_review;
00141 
00142   // different statuses for review
00143   static const unsigned int _remove_heap_entry = 1;
00144   static const unsigned int _review_heap_entry = 2;
00145   static const unsigned int _review_neighbour  = 4;
00146 
00147   /// add a label to a point as to the nature of review needed
00148   /// (includes adding it to list of points needing review) [doesn't
00149   /// affect other labels already set for the point]
00150   void _add_label(Point * point, unsigned int review_flag);
00151 
00152   /// sets the label for the point to be exclusively this 
00153   /// review flag (and adds it to list of points needing review
00154   /// if not already there)
00155   void _set_label(Point * point, unsigned int review_flag);
00156 
00157   /// for all entries of the _points_under_review[] vector, carry out
00158   /// the actions indicated by its review flag; the points are
00159   /// then removed from _points_under_review[] and their flags
00160   /// set to zero
00161   void _deal_with_points_to_review();
00162 
00163   /// carry out the search-tree related operations of point removal
00164   void _remove_from_search_tree(Point * point_to_remove);
00165 
00166   /// carry out the search-tree related operations of point insertion
00167   void _insert_into_search_tree(Point * new_point);
00168 
00169   /// takes a point and creates a shuffle with the given shift
00170   void _point2shuffle(Point & , Shuffle & , unsigned int shift);
00171 
00172   /// pieces needed for converting coordinates to integer
00173   Coord2D _left_corner;
00174   double _range;
00175 
00176   int _ID(const Point *) const;
00177 
00178   triplet<unsigned int> _shifts;     // absolute shifts
00179   triplet<unsigned int> _rel_shifts; // shifts relative to previous shift
00180 
00181   unsigned int _cp_search_range;
00182 };
00183 
00184 
00185 //----------------------------------------------------------------------
00186 /// \if internal_doc
00187 /// @ingroup internal
00188 /// \class ClosestPair2D::Point
00189 /// class for representing all info needed about a point
00190 /// \endif
00191 class ClosestPair2D::Point {
00192 public:
00193   /// the point's coordinates
00194   Coord2D coord;
00195   /// a pointer to its closest neighbour in our structure
00196   Point * neighbour;
00197   /// the corresponding squared distance
00198   double  neighbour_dist2;
00199   /// circulators for each of the shifts of the shuffles
00200   triplet<circulator> circ;
00201 
00202   /// indicates that something special is currently happening to this point
00203   unsigned int review_flag;
00204 
00205   /// returns the distance between two of these objects
00206   double distance2(const Point & other) const {
00207     return coord.distance2(other.coord);
00208   };
00209 
00210   /// creates a shuffle for us with a given shift
00211   //void set_shuffle(Shuffle & shuffle);
00212 };
00213 
00214 
00215 //----------------------------------------------------------------------
00216 /// returns true if floor(ln_base2(x)) < floor(ln_base2(y)), using
00217 /// Chan's neat trick...
00218 inline bool floor_ln2_less(unsigned x, unsigned y) {
00219   if (x>y) return false;
00220   return (x < (x^y)); // beware of operator precedence...
00221 }
00222 
00223 
00224 //----------------------------------------------------------------------
00225 /// returns the ID for the specified point...
00226 inline int ClosestPair2D::_ID(const Point * point) const {
00227   return point - &(_points[0]);
00228 }
00229 
00230 
00231 //
00232 inline unsigned int ClosestPair2D::size() {
00233   return _points.size() - _available_points.size();
00234 }
00235 
00236 
00237 
00238 FASTJET_END_NAMESPACE
00239 
00240 #endif // __FASTJET_CLOSESTPAIR2D__HH__
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends