FastJet  3.0.6
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ClosestPair2D.hh
1 //STARTHEADER
2 // $Id: ClosestPair2D.hh 2577 2011-09-13 15:11:38Z salam $
3 //
4 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5 //
6 //----------------------------------------------------------------------
7 // This file is part of FastJet.
8 //
9 // FastJet is free software; you can redistribute it and/or modify
10 // it under the terms of the GNU General Public License as published by
11 // the Free Software Foundation; either version 2 of the License, or
12 // (at your option) any later version.
13 //
14 // The algorithms that underlie FastJet have required considerable
15 // development and are described in hep-ph/0512210. If you use
16 // FastJet as part of work towards a scientific publication, please
17 // include a citation to the FastJet paper.
18 //
19 // FastJet is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU General Public License
25 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
26 //----------------------------------------------------------------------
27 //ENDHEADER
28 
29 #ifndef __FASTJET_CLOSESTPAIR2D__HH__
30 #define __FASTJET_CLOSESTPAIR2D__HH__
31 
32 #include<vector>
33 #include<stack>
34 #include<iostream>
35 #include "fastjet/internal/ClosestPair2DBase.hh"
36 #include "fastjet/internal/SearchTree.hh"
37 #include "fastjet/internal/MinHeap.hh"
38 
39 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
40 
41 //----------------------------------------------------------------------
42 /// \if internal_doc
43 /// @ingroup internal
44 /// \class ClosestPair2D
45 /// concrete implementation for finding closest pairs in 2D -- will
46 /// use Chan's (hopefully efficient) shuffle based structures
47 /// \endif
48 class ClosestPair2D : public ClosestPair2DBase {
49 public:
50  /// constructor from a vector of 2D positions -- number of objects
51  /// after insertion and deletion must never exceed positions.size();
52  /// objects are given IDs that correspond to their index in the vector
53  /// of positions
54  ClosestPair2D(const std::vector<Coord2D> & positions,
55  const Coord2D & left_corner, const Coord2D & right_corner) {
56  _initialize(positions, left_corner, right_corner, positions.size());
57  };
58 
59  /// constructor which allows structure to grow beyond positions.size(), up
60  /// to max_size
61  ClosestPair2D(const std::vector<Coord2D> & positions,
62  const Coord2D & left_corner, const Coord2D & right_corner,
63  const unsigned int max_size) {
64  _initialize(positions, left_corner, right_corner, max_size);
65  };
66 
67  /// provides the IDs of the closest pair as well as the distance between
68  /// them
69  void closest_pair(unsigned int & ID1, unsigned int & ID2,
70  double & distance2) const;
71 
72  /// removes the entry labelled by ID from the object;
73  void remove(unsigned int ID);
74 
75  /// inserts the position into the closest pair structure and returns the
76  /// ID that has been allocated for the object.
77  unsigned int insert(const Coord2D &);
78 
79  /// removes ID1 and ID2 and inserts position, returning the ID
80  /// corresponding to position...
81  virtual unsigned int replace(unsigned int ID1, unsigned int ID2,
82  const Coord2D & position);
83 
84  /// replaces IDs_to_remove with points at the new_positions
85  /// indicating the IDs allocated to the new points in new_IDs
86  virtual void replace_many(const std::vector<unsigned int> & IDs_to_remove,
87  const std::vector<Coord2D> & new_positions,
88  std::vector<unsigned int> & new_IDs);
89 
90  // mostly for checking how things are working...
91  inline void print_tree_depths(std::ostream & outdev) const {
92  outdev << _trees[0]->max_depth() << " "
93  << _trees[1]->max_depth() << " "
94  << _trees[2]->max_depth() << "\n";
95  };
96 
97  unsigned int size();
98 
99 private:
100 
101  void _initialize(const std::vector<Coord2D> & positions,
102  const Coord2D & left_corner, const Coord2D & right_corner,
103  const unsigned int max_size);
104 
105  static const unsigned int _nshift = 3;
106 
107  class Point; // will be defined below
108 
109  /// since sets of three objects will crop up repeatedly, useful
110  /// to have a triplet class?
111  template<class T> class triplet {
112  public:
113  inline const T & operator[](unsigned int i) const {return _contents[i];};
114  inline T & operator[](unsigned int i) {return _contents[i];};
115  private:
116  T _contents[_nshift];
117  };
118 
119 
120  /// class that will take care of ordering of shuffles for us
121  class Shuffle {
122  public:
123  unsigned int x, y;
124  Point * point;
125  bool operator<(const Shuffle &) const;
126  void operator+=(unsigned int shift) {x += shift; y+= shift;};
127  };
128 
129  typedef SearchTree<Shuffle> Tree;
130  typedef Tree::circulator circulator;
131  typedef Tree::const_circulator const_circulator;
132 
133 
134  triplet<std::auto_ptr<Tree> > _trees;
135  std::auto_ptr<MinHeap> _heap;
136  std::vector<Point> _points;
137  std::stack<Point *> _available_points;
138 
139  /// points that are "under review" in some way
140  std::vector<Point *> _points_under_review;
141 
142  // different statuses for review
143  static const unsigned int _remove_heap_entry = 1;
144  static const unsigned int _review_heap_entry = 2;
145  static const unsigned int _review_neighbour = 4;
146 
147  /// add a label to a point as to the nature of review needed
148  /// (includes adding it to list of points needing review) [doesn't
149  /// affect other labels already set for the point]
150  void _add_label(Point * point, unsigned int review_flag);
151 
152  /// sets the label for the point to be exclusively this
153  /// review flag (and adds it to list of points needing review
154  /// if not already there)
155  void _set_label(Point * point, unsigned int review_flag);
156 
157  /// for all entries of the _points_under_review[] vector, carry out
158  /// the actions indicated by its review flag; the points are
159  /// then removed from _points_under_review[] and their flags
160  /// set to zero
161  void _deal_with_points_to_review();
162 
163  /// carry out the search-tree related operations of point removal
164  void _remove_from_search_tree(Point * point_to_remove);
165 
166  /// carry out the search-tree related operations of point insertion
167  void _insert_into_search_tree(Point * new_point);
168 
169  /// takes a point and creates a shuffle with the given shift
170  void _point2shuffle(Point & , Shuffle & , unsigned int shift);
171 
172  /// pieces needed for converting coordinates to integer
173  Coord2D _left_corner;
174  double _range;
175 
176  int _ID(const Point *) const;
177 
178  triplet<unsigned int> _shifts; // absolute shifts
179  triplet<unsigned int> _rel_shifts; // shifts relative to previous shift
180 
181  unsigned int _cp_search_range;
182 };
183 
184 
185 //----------------------------------------------------------------------
186 /// \if internal_doc
187 /// @ingroup internal
188 /// \class ClosestPair2D::Point
189 /// class for representing all info needed about a point
190 /// \endif
191 class ClosestPair2D::Point {
192 public:
193  /// the point's coordinates
194  Coord2D coord;
195  /// a pointer to its closest neighbour in our structure
196  Point * neighbour;
197  /// the corresponding squared distance
198  double neighbour_dist2;
199  /// circulators for each of the shifts of the shuffles
200  triplet<circulator> circ;
201 
202  /// indicates that something special is currently happening to this point
203  unsigned int review_flag;
204 
205  /// returns the distance between two of these objects
206  double distance2(const Point & other) const {
207  return coord.distance2(other.coord);
208  };
209 
210  /// creates a shuffle for us with a given shift
211  //void set_shuffle(Shuffle & shuffle);
212 };
213 
214 
215 //----------------------------------------------------------------------
216 /// returns true if floor(ln_base2(x)) < floor(ln_base2(y)), using
217 /// Chan's neat trick...
218 inline bool floor_ln2_less(unsigned x, unsigned y) {
219  if (x>y) return false;
220  return (x < (x^y)); // beware of operator precedence...
221 }
222 
223 
224 //----------------------------------------------------------------------
225 /// returns the ID for the specified point...
226 inline int ClosestPair2D::_ID(const Point * point) const {
227  return point - &(_points[0]);
228 }
229 
230 
231 //
232 inline unsigned int ClosestPair2D::size() {
233  return _points.size() - _available_points.size();
234 }
235 
236 
237 
238 FASTJET_END_NAMESPACE
239 
240 #endif // __FASTJET_CLOSESTPAIR2D__HH__
bool floor_ln2_less(unsigned x, unsigned y)
returns true if floor(ln_base2(x)) &lt; floor(ln_base2(y)), using Chan&#39;s neat trick...
bool operator<(SharedPtr< T > const &t, SharedPtr< U > const &u)
comparison: orgering
Definition: SharedPtr.hh:376