FastJet  3.0.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
DnnPlane.cc
1 //STARTHEADER
2 // $Id: DnnPlane.cc 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 
30 #ifndef DROP_CGAL // in case we do not have the code for CGAL
31 
32 #include<set>
33 #include<list>
34 #include "fastjet/internal/DnnPlane.hh"
35 using namespace std;
36 
37 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
38 
39 
40 /// Initialiser from a set of points on an Eta-Phi plane, where both
41 /// eta and phi can have arbitrary ranges
42 DnnPlane::DnnPlane(const vector<EtaPhi> & input_points,
43  const bool & verbose ) {
44 
45  _verbose = verbose;
46  int n = input_points.size();
47 
48  // construct Voronoi diagram in such a way as to get the vertex handles
49  // and remember to set CGAL info with the index of the vertex
50  SuperVertex sv;
51  for (int i = 0; i < n; i++) {
52  sv.vertex =
53  _TR.insert(Point(input_points[i].first, input_points[i].second));
54 
55  // we are not up to dealing with coincident vertices, so make
56  // sure the user knows!
57  _CrashIfVertexPresent(sv.vertex, i);
58 
59  // we need to assicate an index to each vertex -- thus when we get
60  // a vertex (e.g. as a nearest neighbour) from CGAL, we will be
61  // able to figure out which particle it corresponded to.
62  sv.vertex->info() = i;
63  _supervertex.push_back(sv);
64  }
65 
66  // label infinite vertex info with negative index
67  _TR.infinite_vertex()->info() = INFINITE_VERTEX;
68 
69  // set up the structure that holds nearest distances and neighbours
70  for (int j = 0; j < n; j++) {_SetNearest(j);}
71 
72 }
73 
74 
75 //----------------------------------------------------------------------
76 /// Crashes if the given vertex handle already exists. Otherwise
77 /// it does the bookkeeping for future such tests
78 void DnnPlane::_CrashIfVertexPresent(
79  const Vertex_handle & vertex, const int & its_index) {
80  if (!_crash_on_coincidence) return;
81 
82  // vertices that do not have the same geometric position as any
83  // other vertex so far added have info().val() == NEW_VERTEX -- this
84  // is ensured by the InitialisedInt class, which forms the "info"
85  // part of our
86  // CGAL::Triangulation_vertex_base_with_info_2<InitialisedInt,K>
87  //
88  // If the vertex coincides with one that already exists, then
89  // info().val() it's info().val() will have been updated (in
90  // DNN:DNN) to be equal to a vertex "index".
91  if (vertex->info().val() != NEW_VERTEX) {
92  ostringstream err;
93  err << "ERROR in DnnPlane::_CrashIfVertexPresent"
94  <<endl << "Point "<<its_index<<" coincides with point "
95  <<vertex->info().val() << endl;
96  throw DnnError(err.str());
97  }
98 }
99 
100 
101 //----------------------------------------------------------------------
102 /// remove the points labelled by the vector indices_to_remove, and
103 /// add the points specified by the vector points_to_add
104 /// (corresponding indices will be calculated automatically); the
105 /// idea behind this routine is that the points to be added will
106 /// somehow be close to the one or other of the points being removed
107 /// and this can be used by the implementation to provide hints for
108 /// inserting the new points in whatever structure it is using. In a
109 /// kt-algorithm the points being added will be a result of a
110 /// combination of the points to be removed -- hence the proximity
111 /// is (more or less) guaranteed.
112 void DnnPlane::RemoveAndAddPoints(
113  const vector<int> & indices_to_remove,
114  const vector<EtaPhi> & points_to_add,
115  vector<int> & indices_added,
116  vector<int> & indices_of_updated_neighbours) {
117 
118 
119  // build set of UNION of Voronoi neighbours of a pair of nearest
120  // neighbours
121  set<int> NeighbourUnion;
122  // later on it will be convenient to have access to a set (rather
123  // than vector) of indices being removed
124  set<int> indices_removed;
125 
126  // for each of the indices to be removed add the voronoi neighbourhood to
127  // the NeighbourUnion set.
128  for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
129  int index = indices_to_remove[ir];
130  indices_removed.insert(index);
131  if (_verbose) cout << " Starting RemoveAndAddPoints" << endl;
132  if (_verbose) cout << " point " << index << endl;
133  // have a circulators that will go round the Voronoi neighbours of
134  // _supervertex[index1].vertex
135  Vertex_circulator vc = _TR.incident_vertices(_supervertex[index].vertex);
136  Vertex_circulator done = vc;
137  do {
138  // if a neighbouring vertex not the infinite vertex, then add it
139  // to our union of neighbouring vertices.
140  if (_verbose) cout << "examining " << vc->info().val() << endl;
141  if (vc->info().val() != INFINITE_VERTEX) {
142  // NB: from it=1 onwards occasionally it might already have
143  // been inserted -- but double insertion still leaves only one
144  // copy in the set, so there's no problem
145  NeighbourUnion.insert(vc->info().val());
146  if (_verbose) cout << "inserted " << vc->info().val() << endl;
147  }
148  } while (++vc != done);
149  }
150 
151  if (_verbose) {
152  set<int>::iterator it = NeighbourUnion.begin();
153  cout << "Union of neighbours of combined points" << endl;
154  for ( ; it != NeighbourUnion.end(); ++it ) {
155  cout << *it << endl;
156  }
157  }
158 
159  // update set, triangulation and supervertex info
160  for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
161  int index = indices_to_remove[ir];
162 
163  // NeighbourUnion should not contain the points to be removed
164  // (because later we will assume they still exist).
165  NeighbourUnion.erase(indices_to_remove[ir]);
166 
167  // points to be removed should also be eliminated from the
168  // triangulation and the supervertex structure should be updated
169  // to reflect the fact that the points are no longer valid.
170  _TR.remove(_supervertex[index].vertex);
171  _supervertex[index].vertex = NULL;
172  }
173 
174  // add new point: give a "hint" to the inserter that
175  // the new point should be added close to old points -- the easiest way
176  // of getting this is to take a point from the NeighbourUnion AFTER we have
177  // removed point1, point2, and to get one of its incident faces.
178  //
179  // This hinting improves speed by c. 25% for 10^4 points because it
180  // avoids the need for a costly (sqrt{N}) location step (at least
181  // with a non-hierarchical triangulation -- with a hierarchical one,
182  // this step could be done away with, though there will still be a
183  // cost of O(ln N) to pay.
184  //
185  // For some reason inserting the point before the two removals
186  // slows things down by c. 25%. This importance of the order
187  // is not understood.
188  //
189  // At some point it might be worth trying to select the "nearest"
190  // of the various points in the neighbour union to avoid large
191  // steps in cases where we have 0..2pi periodicity and the first member
192  // of the neighbour union happens to be on the wrong side.
193  Face_handle face;
194  if (indices_to_remove.size() > 0) {
195  // face can only be found if there were points to remove in first place
196  face = _TR.incident_faces(
197  _supervertex[*NeighbourUnion.begin()].vertex);}
198  // make sure the output arrays are empty
199  indices_added.clear();
200  indices_of_updated_neighbours.clear();
201  for (size_t ia = 0; ia < points_to_add.size(); ia++) {
202  SuperVertex sv;
203  _supervertex.push_back(sv);
204  int index = _supervertex.size()-1;
205  indices_added.push_back(index);
206 
207  if (indices_to_remove.size() > 0) {
208  // be careful of using face (for location hinting) only when it exists
209  _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first,
210  points_to_add[ia].second),face);}
211  else {
212  _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first,
213  points_to_add[ia].second));
214  }
215  // we are not up to dealing with coincident vertices, so make
216  // sure the user knows!
217  _CrashIfVertexPresent(_supervertex[index].vertex, index);
218  _supervertex[index].vertex->info() = index;
219 
220  // first find nearest neighbour of "newpoint" (shorthand for
221  // _supervertex[index].vertex); while we're at it, for each of the
222  // voronoi neighbours, "D", of newpoint, examine whether newpoint is
223  // closer to "D" than D's current nearest neighbour -- when this
224  // occurs, put D into indices_of_updated_neighbours.
225  //
226  // manually put newpoint on indices_of_updated_neighbours
227  indices_of_updated_neighbours.push_back(index);
228  _SetAndUpdateNearest(index, indices_of_updated_neighbours);
229  }
230 
231  // for Voronoi neighbours j of any of the removed points for which
232  // one of those removed points was the nearest neighbour,
233  // redetermine the nearest neighbour of j and add j onto the vector
234  // of indices_of_updated_neighbours.
235  set<int>::iterator it2 = NeighbourUnion.begin();
236  for ( ; it2 != NeighbourUnion.end(); ++it2 ) {
237  int j = *it2;
238  // the if avoids the vertex at infinity, which gets a negative index
239  if( j != INFINITE_VERTEX ) {
240  // this is where we check if the nearest neighbour of j was one
241  // of the removed points
242  if (indices_removed.count(_supervertex[j].NNindex)) {
243  if (_verbose) cout << "j " << j << endl;
244  _SetNearest(j);
245  indices_of_updated_neighbours.push_back(j);
246  if (_verbose) cout << "NN of " << j << " : "
247  << _supervertex[j].NNindex
248  << ", dist = " << _supervertex[j].NNdistance <<endl;
249  }
250  }
251  }
252 
253 }
254 
255 
256 //----------------------------------------------------------------------
257 /// Determines the index and distance of the nearest neighbour to
258 /// point j and puts the information into the _supervertex entry for j.
259 void DnnPlane::_SetNearest (const int & j) {
260  Vertex_handle current = _supervertex[j].vertex;
261  Vertex_circulator vc = _TR.incident_vertices(current);
262  Vertex_circulator done = vc;
263  double dist;
264  double mindist = HUGE_DOUBLE; // change this to "HUGE" or max_double?
265  Vertex_handle nearest = _TR.infinite_vertex();
266 
267  // when there is only one finite point left in the triangulation,
268  // there are no triangles. Presumably this is why voronoi returns
269  // NULL for the incident vertex circulator. Check if this is
270  // happening before circulating over it... (Otherwise it crashes
271  // when looking for neighbours of last point)
272  if (vc != NULL) do {
273  if ( vc->info().val() != INFINITE_VERTEX) {
274  // find distance between j and its Voronoi neighbour (vc)
275  if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
276  dist = _euclid_distance(current->point(), vc->point());
277  // check if j is closer to vc than vc's currently registered
278  // nearest neighbour (and update things if it is)
279  if (dist < mindist) {
280  mindist = dist; nearest = vc;
281  if (_verbose) cout << "nearer ";
282  }
283  if (_verbose) cout << vc->point() << "; "<< dist << endl;
284  }
285  } while (++vc != done); // move on to next Voronoi neighbour
286  // set j's supervertex info about nearest neighbour
287  _supervertex[j].NNindex = nearest->info().val();
288  _supervertex[j].NNdistance = mindist;
289 }
290 
291 //----------------------------------------------------------------------
292 /// Determines and stores the nearest neighbour of j, and where
293 /// necessary update the nearest-neighbour info of Voronoi neighbours
294 /// of j;
295 ///
296 /// For each voronoi neighbour D of j if the distance between j and D
297 /// is less than D's own nearest neighbour, then update the
298 /// nearest-neighbour info in D; push D's index onto
299 /// indices_of_updated_neighbours
300 ///
301 /// Note that j is NOT pushed onto indices_of_updated_neighbours --
302 /// if you want it there, put it there yourself.
303 ///
304 /// NB: note that we have _SetAndUpdateNearest as a completely
305 /// separate routine from _SetNearest because we want to
306 /// use one single ciruclation over voronoi neighbours to find the
307 /// nearest neighbour and to update the voronoi neighbours if need
308 /// be.
309 void DnnPlane::_SetAndUpdateNearest(
310  const int & j,
311  vector<int> & indices_of_updated_neighbours) {
312 
313  Vertex_handle current = _supervertex[j].vertex;
314  Vertex_circulator vc = _TR.incident_vertices(current);
315  Vertex_circulator done = vc;
316  double dist;
317  double mindist = HUGE_DOUBLE; // change this to "HUGE" or max_double?
318  Vertex_handle nearest = _TR.infinite_vertex();
319 
320  // when there is only one finite point left in the triangulation,
321  // there are no triangles. Presumably this is why voronoi returns
322  // NULL for the incident vertex circulator. Check if this is
323  // happening before circulating over it... (Otherwise it crashes
324  // when looking for neighbours of last point)
325  if (vc != NULL) do {
326  if (vc->info().val() != INFINITE_VERTEX) {
327  if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
328  // find distance between j and its Voronoi neighbour (vc)
329  dist = _euclid_distance(current->point(), vc->point());
330  // update the mindist if we are closer than anything found so far
331  if (dist < mindist) {
332  mindist = dist; nearest = vc;
333  if (_verbose) cout << "nearer ";
334  }
335  // find index corresponding to vc for easy manipulation
336  int vcindx = vc->info().val();
337  if (_verbose) cout << vc->point() << "; "<< dist << endl;
338  // check if j is closer to vc than vc's currently registered
339  // nearest neighbour (and update things if it is)
340  if (dist < _supervertex[vcindx].NNdistance) {
341  if (_verbose) cout << vcindx << "'s NN becomes " << current->info().val() << endl;
342  _supervertex[vcindx].NNdistance = dist;
343  _supervertex[vcindx].NNindex = j;
344  indices_of_updated_neighbours.push_back(vcindx);
345  }
346  }
347  } while (++vc != done); // move on to next Voronoi neighbour
348  // set j's supervertex info about nearest neighbour
349  _supervertex[j].NNindex = nearest->info().val();
350  _supervertex[j].NNdistance = mindist;
351 }
352 
353 FASTJET_END_NAMESPACE
354 
355 #endif // DROP_CGAL