FastJet  3.1.2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
DnnPlane.cc
1 //FJSTARTHEADER
2 // $Id: DnnPlane.cc 3442 2014-07-24 07:20:49Z salam $
3 //
4 // Copyright (c) 2005-2014, 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 
32 #ifndef DROP_CGAL // in case we do not have the code for CGAL
33 
34 #include<set>
35 #include<list>
36 #include "fastjet/internal/DnnPlane.hh"
37 
38 using namespace std;
39 
40 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
41 
42 const double DnnPlane::DISTANCE_FOR_CGAL_CHECKS=1.0e-12;
43 
44 
45 /// Initialiser from a set of points on an Eta-Phi plane, where both
46 /// eta and phi can have arbitrary ranges
47 DnnPlane::DnnPlane(const vector<EtaPhi> & input_points,
48  const bool & verbose ) {
49 
50  _verbose = verbose;
51  int n = input_points.size();
52 
53  // construct Voronoi diagram in such a way as to get the vertex handles
54  // and remember to set CGAL info with the index of the vertex
55  SuperVertex sv;
56  for (int i = 0; i < n; i++) {
57  sv.vertex =
58  _TR.insert(Point(input_points[i].first, input_points[i].second));
59 
60  // check if we are dealing with coincident vertices
61  int coinciding_index = _CheckIfVertexPresent(sv.vertex, i);
62  if (coinciding_index == i){
63  // we need to associate an index to each vertex -- thus when we get
64  // a vertex (e.g. as a nearest neighbour) from CGAL, we will be
65  // able to figure out which particle it corresponded to.
66  sv.vertex->info() = sv.coincidence = i;
67  } else {
68  //cout << " coincident with " << coinciding_index << endl;
69  // the new vertex points to the already existing one and we
70  // record the coincidence
71  //
72  // Note that we must not only set the coincidence of the
73  // currently-added particle, the one it coincides with also
74  // needs be updated (taking into account that it might already
75  // coincide with another one)
76  //
77  // An example may help. Say coinciding_index = i1 and we're adding i2==i.
78  // Then _sv[i2].coincidence = i1; _sv[i1].coincidence = i2. In both
79  // cases sv.vertex->info() == i1;
80  //
81  // Later on we add i3; we find out that its coinciding index is i1;
82  // so we set _sv[i3].coincidence = i2 and sv[i1].coincidence = i3.
83  //
84  // This gives us the structure
85  // _supervertex[i1].coincidence == in
86  // _supervertex[i2].coincidence == i1
87  // ...
88  // _supervertex[in].coincidence == in-1
89  //
90  sv.coincidence = _supervertex[coinciding_index].coincidence; // handles cases with previous coincidences
91  _supervertex[coinciding_index].coincidence = i;
92  }
93 
94  _supervertex.push_back(sv);
95  }
96 
97  // label infinite vertex info with negative index
98  _TR.infinite_vertex()->info() = INFINITE_VERTEX;
99 
100  // set up the structure that holds nearest distances and neighbours
101  for (int j = 0; j < n; j++) {_SetNearest(j);}
102 
103 }
104 
105 
106 //----------------------------------------------------------------------
107 /// Crashes if the given vertex handle already exists. Otherwise
108 /// it does the bookkeeping for future such tests
109 int DnnPlane::_CheckIfVertexPresent(
110  const Vertex_handle & vertex, const int its_index) {
111  // vertices that do not have the same geometric position as any
112  // other vertex so far added have info().val() == NEW_VERTEX -- this
113  // is ensured by the InitialisedInt class, which forms the "info"
114  // part of our
115  // CGAL::Triangulation_vertex_base_with_info_2<InitialisedInt,K>
116  //
117  // If the vertex coincides with one that already exists, then
118  // info().val() it's info().val() will have been updated (in
119  // DNN:DNN) to be equal to a vertex "index".
120  if (vertex->info().val() != NEW_VERTEX) {
121  if (_crash_on_coincidence){
122  ostringstream err;
123  err << "Error: DnnPlane::_CheckIfVertexPresent"
124  << "Point "<<its_index<<" coincides with point "
125  <<vertex->info().val() << endl;
126  throw DnnError(err.str());
127  }
128  return vertex->info().val();
129  }
130 
131  return its_index;
132 }
133 
134 
135 //----------------------------------------------------------------------
136 /// remove the points labelled by the vector indices_to_remove, and
137 /// add the points specified by the vector points_to_add
138 /// (corresponding indices will be calculated automatically); the
139 /// idea behind this routine is that the points to be added will
140 /// somehow be close to the one or other of the points being removed
141 /// and this can be used by the implementation to provide hints for
142 /// inserting the new points in whatever structure it is using. In a
143 /// kt-algorithm the points being added will be a result of a
144 /// combination of the points to be removed -- hence the proximity
145 /// is (more or less) guaranteed.
146 void DnnPlane::RemoveAndAddPoints(
147  const vector<int> & indices_to_remove,
148  const vector<EtaPhi> & points_to_add,
149  vector<int> & indices_added,
150  vector<int> & indices_of_updated_neighbours) {
151 
152  if (_verbose) cout << "Starting DnnPlane::RemoveAndAddPoints" << endl;
153 
154  // build set of UNION of Voronoi neighbours of a pair of nearest
155  // neighbours
156  set<int> NeighbourUnion;
157  // later on it will be convenient to have access to a set (rather
158  // than vector) of indices being removed
159  set<int> indices_removed;
160 
161  // for each of the indices to be removed add the voronoi
162  // neighbourhood to the NeighbourUnion set as well as the coinciding
163  // points that had the current point as coincidence before.
164  for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
165  int index = indices_to_remove[ir];
166  indices_removed.insert(index);
167  if (_verbose) cout << " scheduling point " << index << " for removal" << endl;
168 
169  if (_supervertex[index].coincidence != index){
170  // we have a coincidence
171  //
172  // The only one of the coincident points that has to be
173  // inserted in the neighbourhood list (and thus updated) is the
174  // one that has 'index' as coincidence.
175  int new_index = _supervertex[index].coincidence;
176  while (_supervertex[new_index].coincidence != index)
177  new_index = _supervertex[new_index].coincidence;
178  if (_verbose) cout << " inserted coinciding " << new_index << " to neighbours union" << endl;
179  NeighbourUnion.insert(new_index);
180 
181  // if this is the point among the coiciding ones that holds the
182  // CGAL vertex, then also insert the CGAL neighbours, otherwise
183  // just skip that step.
184  if (index != _supervertex[index].vertex->info().val()) continue;
185  }
186 
187  // have a circulators that will go round the Voronoi neighbours of
188  // _supervertex[index1].vertex
189  Vertex_circulator vc = _TR.incident_vertices(_supervertex[index].vertex);
190  Vertex_circulator done = vc;
191  if (vc != NULL){ // a safety check in case there is no Voronoi
192  // neighbour (which may happen e.g. if we just
193  // have a bunch of coincident points)
194  do {
195  // if a neighbouring vertex is not the infinite vertex, then add it
196  // to our union of neighbouring vertices.
197  if (_verbose) cout << "examining " << vc->info().val() << endl;
198  if (vc->info().val() != INFINITE_VERTEX) {
199  // NB: from it=1 onwards occasionally it might already have
200  // been inserted -- but double insertion still leaves only one
201  // copy in the set, so there's no problem
202  NeighbourUnion.insert(vc->info().val());
203  if (_verbose) cout << " inserted " << vc->info().val() << " to neighbours union" << endl;
204  }
205  } while (++vc != done);
206  }
207  }
208 
209  if (_verbose) {
210  set<int>::iterator it = NeighbourUnion.begin();
211  cout << "Union of neighbours of combined points" << endl;
212  for ( ; it != NeighbourUnion.end(); ++it ) {
213  cout << *it << endl;
214  }
215  }
216 
217  // update set, triangulation and supervertex info
218  for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
219  int index = indices_to_remove[ir];
220  if (_verbose) cout << " removing " << index << endl;
221 
222  // NeighbourUnion should not contain the points to be removed
223  // (because later we will assume they still exist).
224  NeighbourUnion.erase(indices_to_remove[ir]);
225 
226  // first deal with coincidences
227  if (_supervertex[index].coincidence != index){
228  int new_index = _supervertex[index].coincidence;
229 
230  // if this is the point among the coiciding ones that "owns" the
231  // CGAL vertex we need to re-label the CGAL vertex so that it
232  // points to the coincident particle and set the current one to
233  // NULL
234  //
235  // This can be done only on the first point as they all share
236  // the same value
237  //
238  // Note that this has to be done before the following step since
239  // it will alter the coincidence information
240  if (index == _supervertex[index].vertex->info().val())
241  _supervertex[new_index].vertex->info() = new_index;
242 
243  // we need to browse the coincidences until we end the loop, at
244  // which point we reset the coincidence of the point that has
245  // the current one as a coincidence
246  while (_supervertex[new_index].coincidence != index)
247  new_index = _supervertex[new_index].coincidence;
248  _supervertex[new_index].coincidence = _supervertex[index].coincidence;
249 
250  // remove the coincidence on the point being removed and mark it
251  // as removed
252  _supervertex[index].coincidence = index;
253  _supervertex[index].vertex = NULL;
254 
255  continue;
256  }
257 
258  // points to be removed should also be eliminated from the
259  // triangulation and the supervertex structure should be updated
260  // to reflect the fact that the points are no longer valid.
261  _TR.remove(_supervertex[index].vertex);
262  _supervertex[index].vertex = NULL;
263  }
264 
265  // add new point: give a "hint" to the inserter that
266  // the new point should be added close to old points -- the easiest way
267  // of getting this is to take a point from the NeighbourUnion AFTER we have
268  // removed point1, point2, and to get one of its incident faces.
269  //
270  // This hinting improves speed by c. 25% for 10^4 points because it
271  // avoids the need for a costly (sqrt{N}) location step (at least
272  // with a non-hierarchical triangulation -- with a hierarchical one,
273  // this step could be done away with, though there will still be a
274  // cost of O(ln N) to pay.
275  //
276  // For some reason inserting the point before the two removals
277  // slows things down by c. 25%. This importance of the order
278  // is not understood.
279  //
280  // At some point it might be worth trying to select the "nearest"
281  // of the various points in the neighbour union to avoid large
282  // steps in cases where we have 0..2pi periodicity and the first member
283  // of the neighbour union happens to be on the wrong side.
284  Face_handle face;
285  //if (indices_to_remove.size() > 0) { // GS: use NeighbourUnion instead
286  // (safe also in case of coincidences)
287  if (NeighbourUnion.size() > 0) {
288  // face can only be found if there were points to remove in first place
289  face = _TR.incident_faces(
290  _supervertex[*NeighbourUnion.begin()].vertex);}
291  // make sure the output arrays are empty
292  indices_added.clear();
293  indices_of_updated_neighbours.clear();
294  for (size_t ia = 0; ia < points_to_add.size(); ia++) {
295  SuperVertex sv;
296  _supervertex.push_back(sv);
297  int index = _supervertex.size()-1;
298  indices_added.push_back(index);
299  if (_verbose) cout << " adding " << index << endl;
300 
301  //if (indices_to_remove.size() > 0) {
302  if (NeighbourUnion.size() > 0) {
303  // be careful of using face (for location hinting) only when it exists
304  _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first,
305  points_to_add[ia].second),face);}
306  else {
307  _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first,
308  points_to_add[ia].second));
309  }
310 
311  // check if this leads to a coincidence
312  int coinciding_index = _CheckIfVertexPresent(_supervertex[index].vertex, index);
313  if (coinciding_index == index){
314  // we need to associate an index to each vertex -- thus when we get
315  // a vertex (e.g. as a nearest neighbour) from CGAL, we will be
316  // able to figure out which particle it corresponded to.
317  _supervertex[index].vertex->info() = _supervertex[index].coincidence = index;
318  } else {
319  if (_verbose) cout << " coinciding with vertex " << coinciding_index << endl;
320  // the new vertex points to an already existing one and we
321  // record the coincidence
322  //
323  // we also update the NN of the coinciding particle (to avoid
324  // having to loop over the list of coinciding neighbours later)
325  // This is done first as it allows us to check if this is a new
326  // coincidence or a coincidence added to a particle that was
327  // previously "alone"
328  _supervertex[coinciding_index].NNindex = index;
329  _supervertex[coinciding_index].NNdistance = 0.0;
330  indices_of_updated_neighbours.push_back(coinciding_index);
331 
332  // Note that we must not only set the coincidence of the
333  // currently-added particle, the one it coincides with also
334  // needs be updated (taking into account that it might already
335  // coincide with another one)
336  _supervertex[index].coincidence = _supervertex[coinciding_index].coincidence; // handles cases with previous coincidences
337  _supervertex[coinciding_index].coincidence = index;
338 
339  }
340 
341  // first find nearest neighbour of "newpoint" (shorthand for
342  // _supervertex[index].vertex); while we're at it, for each of the
343  // voronoi neighbours, "D", of newpoint, examine whether newpoint is
344  // closer to "D" than D's current nearest neighbour -- when this
345  // occurs, put D into indices_of_updated_neighbours.
346  //
347  // manually put newpoint on indices_of_updated_neighbours
348  indices_of_updated_neighbours.push_back(index);
349  _SetAndUpdateNearest(index, indices_of_updated_neighbours);
350 
351  //cout << "Added: " << setprecision(20) << " ("
352  // << points_to_add[ia].first << "," << points_to_add[ia].second
353  // << ") with index " << index << endl;
354  }
355 
356  // for Voronoi neighbours j of any of the removed points for which
357  // one of those removed points was the nearest neighbour,
358  // redetermine the nearest neighbour of j and add j onto the vector
359  // of indices_of_updated_neighbours.
360  set<int>::iterator it2 = NeighbourUnion.begin();
361  for ( ; it2 != NeighbourUnion.end(); ++it2 ) {
362  int j = *it2;
363  // the if avoids the vertex at infinity, which gets a negative index
364  if( j != INFINITE_VERTEX ) {
365  // this is where we check if the nearest neighbour of j was one
366  // of the removed points
367  if (indices_removed.count(_supervertex[j].NNindex)) {
368  if (_verbose) cout << "j " << j << endl;
369  _SetNearest(j);
370  indices_of_updated_neighbours.push_back(j);
371  if (_verbose) cout << "NN of " << j << " : "
372  << _supervertex[j].NNindex
373  << ", dist = " << _supervertex[j].NNdistance <<endl;
374  }
375  }
376  }
377 
378  if (_verbose) cout << "Leaving DnnPlane::RemoveAndAddPoints" << endl;
379 }
380 
381 //----------------------------------------------------------------------
382 /// Determines the index and distance of the nearest neighbour to
383 /// point j and puts the information into the _supervertex entry for j.
384 void DnnPlane::_SetNearest (const int j) {
385  // first deal with the cases where we have a coincidence
386  if (_supervertex[j].coincidence != j){
387  _supervertex[j].NNindex = _supervertex[j].coincidence;
388  _supervertex[j].NNdistance = 0.0;
389  return;
390  }
391 
392  // The code below entirely uses CGAL distance comparisons to compute
393  // the nearest neighbour. It has the mais drawback to induice a
394  // 10-20% time penalty so we switched to our own comparison (which
395  // only turns to CGAL for dangerous situations)
396  //
397  // Vertex_handle current = _supervertex[j].vertex;
398  // Vertex_circulator vc = _TR.incident_vertices(current);
399  // Vertex_circulator done = vc;
400  // Vertex_handle nearest = _TR.infinite_vertex();
401  // double mindist = HUGE_DOUBLE;
402  //
403  // // when there is only one finite point left in the triangulation,
404  // // there are no triangles. Presumably this is why voronoi returns
405  // // NULL for the incident vertex circulator. Check if this is
406  // // happening before circulating over it... (Otherwise it crashes
407  // // when looking for neighbours of last point)
408  // if (vc != NULL){
409  // // initialise the nearest vertex handle to the first incident
410  // // vertex that is not INFINITE_VERTEX
411  // while (vc->info().val() == INFINITE_VERTEX){
412  // vc++;
413  // if (vc==done) break; // if vc==done, then INFINITE_VERTEX is the
414  // // only element in the neighbourhood
415  // }
416  //
417  // // if there is just the infinite vertex, we have vc->info().val()
418  // // == INFINITE_VERTEX and nothing has to be done
419  // // otherwise, use the current vc as an initialisation
420  // if (vc->info().val() != INFINITE_VERTEX){
421  // nearest = vc; // initialisation to the first non-infinite vertex
422  //
423  // // and loop over the following ones
424  // while (++vc != done){
425  // // we should not compare with the infinite vertex
426  // if (vc->info().val() == INFINITE_VERTEX) continue;
427  //
428  // if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
429  // // use CGAL's distance comparison to check if 'vc' is closer to
430  // // 'current' than the nearest so far (we include the == case for
431  // // safety though it should not matter in this precise case)
432  // if (CGAL::compare_distance_to_point(current->point(), vc->point(), nearest->point())!=CGAL::LARGER){
433  // nearest = vc;
434  // if (_verbose) cout << "nearer";
435  // }
436  // }
437  //
438  // // now compute the distance
439  // //
440  // // Note that since we're always using CGAL to compare distances
441  // // (and never the distance computed using _euclid_distance) we
442  // // should not worry about rounding errors in mindist
443  // mindist = _euclid_distance(current->point(), nearest->point());
444  // }
445  // }
446  //
447  // // set j's supervertex info about nearest neighbour
448  // _supervertex[j].NNindex = nearest->info().val();
449  // _supervertex[j].NNdistance = mindist;
450 
451  Vertex_handle current = _supervertex[j].vertex;
452  Vertex_circulator vc = _TR.incident_vertices(current);
453  Vertex_circulator done = vc;
454  double dist;
455  double mindist = HUGE_DOUBLE; // change this to "HUGE" or max_double?
456  Vertex_handle nearest = _TR.infinite_vertex();
457 
458  // when there is only one finite point left in the triangulation,
459  // there are no triangles. Presumably this is why voronoi returns
460  // NULL for the incident vertex circulator. Check if this is
461  // happening before circulating over it... (Otherwise it crashes
462  // when looking for neighbours of last point)
463  if (vc != NULL) do {
464  if ( vc->info().val() != INFINITE_VERTEX) {
465  // find distance between j and its Voronoi neighbour (vc)
466  if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
467 
468  // check if j is closer to vc than vc's currently registered
469  // nearest neighbour (and update things if it is)
470  if (_is_closer_to(current->point(), vc->point(), nearest, dist, mindist)){
471  nearest = vc;
472  if (_verbose) cout << "nearer ";
473  }
474  if (_verbose) cout << vc->point() << "; "<< dist << endl;
475  }
476  } while (++vc != done); // move on to next Voronoi neighbour
477 
478  // set j's supervertex info about nearest neighbour
479  _supervertex[j].NNindex = nearest->info().val();
480  _supervertex[j].NNdistance = mindist;
481 }
482 
483 //----------------------------------------------------------------------
484 /// Determines and stores the nearest neighbour of j, and where
485 /// necessary updates the nearest-neighbour info of Voronoi neighbours
486 /// of j;
487 ///
488 /// For each voronoi neighbour D of j if the distance between j and D
489 /// is less than D's own nearest neighbour, then update the
490 /// nearest-neighbour info in D; push D's index onto
491 /// indices_of_updated_neighbours
492 ///
493 /// Note that j is NOT pushed onto indices_of_updated_neighbours --
494 /// if you want it there, put it there yourself.
495 ///
496 /// NB: note that we have _SetAndUpdateNearest as a completely
497 /// separate routine from _SetNearest because we want to
498 /// use one single circulation over voronoi neighbours to find the
499 /// nearest neighbour and to update the voronoi neighbours if need
500 /// be.
501 void DnnPlane::_SetAndUpdateNearest(
502  const int j,
503  vector<int> & indices_of_updated_neighbours) {
504  //cout << "SetAndUpdateNearest for point " << j << endl;
505  // first deal with coincidences
506  if (_supervertex[j].coincidence != j){
507  _supervertex[j].NNindex = _supervertex[j].coincidence;
508  _supervertex[j].NNdistance = 0.0;
509  //cout << " set to coinciding point " << _supervertex[j].coincidence << endl;
510  return;
511  }
512 
513  Vertex_handle current = _supervertex[j].vertex;
514  Vertex_circulator vc = _TR.incident_vertices(current);
515  Vertex_circulator done = vc;
516  double dist;
517  double mindist = HUGE_DOUBLE; // change this to "HUGE" or max_double?
518  Vertex_handle nearest = _TR.infinite_vertex();
519 
520  // when there is only one finite point left in the triangulation,
521  // there are no triangles. Presumably this is why voronoi returns
522  // NULL for the incident vertex circulator. Check if this is
523  // happening before circulating over it... (Otherwise it crashes
524  // when looking for neighbours of last point)
525  if (vc != NULL) do {
526  if (vc->info().val() != INFINITE_VERTEX) {
527  if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
528 
529  // update the mindist if we are closer than anything found so far
530  if (_is_closer_to(current->point(), vc->point(), nearest, dist, mindist)){
531  nearest = vc;
532  if (_verbose) cout << "nearer ";
533  }
534 
535  // find index corresponding to vc for easy manipulation
536  int vcindx = vc->info().val();
537  if (_verbose) cout << vc->point() << "; "<< dist << endl;
538 
539  if (_is_closer_to_with_hint(vc->point(), current->point(),
540  _supervertex[_supervertex[vcindx].NNindex].vertex,
541  dist, _supervertex[vcindx].NNdistance)){
542  if (_verbose) cout << vcindx << "'s NN becomes " << current->info().val() << endl;
543  _supervertex[vcindx].NNindex = j;
544  indices_of_updated_neighbours.push_back(vcindx);
545  }
546 
547  // original code without the use of CGAL distance in potentially
548  // dangerous cases
549  //
550  // // check if j is closer to vc than vc's currently registered
551  // // nearest neighbour (and update things if it is)
552  // //
553  // // GS: originally, the distance test below was a strict <. It
554  // // has to be <= because if the two distances are ==, it is
555  // // possible that the old NN is no longer connected to vc in
556  // // the triangulation, and we are sure that the newly
557  // // inserted point (j) is (since we loop over j's
558  // // neighbouring points in the triangulation).
559  // if (dist <= _supervertex[vcindx].NNdistance) {
560  // if (_verbose) cout << vcindx << "'s NN becomes " << current->info().val() << endl;
561  // _supervertex[vcindx].NNdistance = dist;
562  // _supervertex[vcindx].NNindex = j;
563  // indices_of_updated_neighbours.push_back(vcindx);
564  // }
565  }
566  } while (++vc != done); // move on to next Voronoi neighbour
567  // set j's supervertex info about nearest neighbour
568  //cout << " set to point " << nearest->info().val() << endl;
569  _supervertex[j].NNindex = nearest->info().val();
570  _supervertex[j].NNdistance = mindist;
571 }
572 
573 FASTJET_END_NAMESPACE
574 
575 #endif // DROP_CGAL
Triangulation::Vertex_circulator Vertex_circulator
CGAL Point structure.