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