FastJet 3.0.2
|
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__