Index: external/fastjet/DnnPlane.cc
===================================================================
--- external/fastjet/DnnPlane.cc	(revision d7d2da362406b9e37fdf238970417c94779b54b0)
+++ external/fastjet/DnnPlane.cc	(revision be2222cc14ff474f18d30a4e0ebc21a952777df5)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id$
+//FJSTARTHEADER
+// $Id: DnnPlane.cc 3442 2014-07-24 07:20:49Z salam $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  development. They are described in the original FastJet paper,
+//  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
 //  FastJet as part of work towards a scientific publication, please
-//  include a citation to the FastJet paper.
+//  quote the version you use and include a citation to the manual and
+//  optionally also to hep-ph/0512210.
 //
 //  FastJet is distributed in the hope that it will be useful,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 
@@ -33,7 +35,10 @@
 #include<list>
 #include "fastjet/internal/DnnPlane.hh"
+
 using namespace std;
 
 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+const double DnnPlane::DISTANCE_FOR_CGAL_CHECKS=1.0e-12;  
 
 
@@ -53,12 +58,38 @@
        _TR.insert(Point(input_points[i].first, input_points[i].second));
 
-    // we are not up to dealing with coincident vertices, so make 
-    // sure the user knows!
-    _CrashIfVertexPresent(sv.vertex, i);
-    
-    // we need to assicate an index to each vertex -- thus when we get
-    // a vertex (e.g. as a nearest neighbour) from CGAL, we will be
-    // able to figure out which particle it corresponded to.
-    sv.vertex->info() = i;
+    // check if we are dealing with coincident vertices
+    int coinciding_index = _CheckIfVertexPresent(sv.vertex, i);
+    if (coinciding_index == i){
+      // we need to associate an index to each vertex -- thus when we get
+      // a vertex (e.g. as a nearest neighbour) from CGAL, we will be
+      // able to figure out which particle it corresponded to.
+      sv.vertex->info() = sv.coincidence = i;
+    } else {
+      //cout << "  coincident with " << coinciding_index << endl;
+      // the new vertex points to the already existing one and we
+      // record the coincidence
+      //
+      // Note that we must not only set the coincidence of the
+      // currently-added particle, the one it coincides with also
+      // needs be updated (taking into account that it might already
+      // coincide with another one)
+      //
+      // An example may help. Say coinciding_index = i1 and we're adding i2==i.
+      // Then _sv[i2].coincidence = i1; _sv[i1].coincidence = i2. In both
+      // cases sv.vertex->info() == i1;
+      //
+      // Later on we add i3; we find out that its coinciding index is i1;
+      // so we set _sv[i3].coincidence = i2 and sv[i1].coincidence = i3. 
+      //
+      // This gives us the structure 
+      //  _supervertex[i1].coincidence == in
+      //  _supervertex[i2].coincidence == i1
+      //  ...
+      //  _supervertex[in].coincidence == in-1
+      //
+      sv.coincidence = _supervertex[coinciding_index].coincidence; // handles cases with previous coincidences
+      _supervertex[coinciding_index].coincidence = i;
+    }
+      
     _supervertex.push_back(sv);    
   }
@@ -76,8 +107,6 @@
 /// Crashes if the given vertex handle already exists. Otherwise
 /// it does the bookkeeping for future such tests
-void DnnPlane::_CrashIfVertexPresent(
-	const Vertex_handle & vertex, const int & its_index) {
-  if (!_crash_on_coincidence) return;
-
+int DnnPlane::_CheckIfVertexPresent(
+	const Vertex_handle & vertex, const int its_index) {
   // vertices that do not have the same geometric position as any
   // other vertex so far added have info().val() == NEW_VERTEX -- this
@@ -90,10 +119,15 @@
   // DNN:DNN) to be equal to a vertex "index".
   if (vertex->info().val() != NEW_VERTEX) {
-    ostringstream err;
-    err << "ERROR in DnnPlane::_CrashIfVertexPresent"
-	 <<endl << "Point "<<its_index<<" coincides with point "
-	 <<vertex->info().val() << endl;
-    throw DnnError(err.str());
-  } 
+    if (_crash_on_coincidence){
+      ostringstream err;
+      err << "Error: DnnPlane::_CheckIfVertexPresent"
+	  << "Point "<<its_index<<" coincides with point "
+	  <<vertex->info().val() << endl;
+      throw DnnError(err.str());
+    }
+    return vertex->info().val();
+  }
+
+  return its_index;
 }
 
@@ -116,4 +150,5 @@
 			  vector<int> & indices_of_updated_neighbours) {
 
+  if (_verbose) cout << "Starting  DnnPlane::RemoveAndAddPoints" << endl;
 
   // build set of UNION of Voronoi neighbours of a pair of nearest
@@ -124,27 +159,50 @@
   set<int> indices_removed;
 
-  // for each of the indices to be removed add the voronoi neighbourhood to
-  // the NeighbourUnion set.
+  // for each of the indices to be removed add the voronoi
+  // neighbourhood to the NeighbourUnion set as well as the coinciding
+  // points that had the current point as coincidence before.
   for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
     int index = indices_to_remove[ir];
     indices_removed.insert(index);
-    if (_verbose) cout << " Starting  RemoveAndAddPoints" << endl; 
-    if (_verbose) cout << " point " << index << endl;			 
+    if (_verbose) cout << "  scheduling point " << index << " for removal" << endl;
+
+    if (_supervertex[index].coincidence != index){
+      // we have a coincidence
+      //
+      // The only one of the coincident points that has to be
+      // inserted in the neighbourhood list (and thus updated) is the
+      // one that has 'index' as coincidence.
+      int new_index = _supervertex[index].coincidence;
+      while (_supervertex[new_index].coincidence != index)
+	new_index = _supervertex[new_index].coincidence;
+      if (_verbose) cout << "  inserted coinciding " << new_index << " to neighbours union" << endl;
+      NeighbourUnion.insert(new_index);
+
+      // if this is the point among the coiciding ones that holds the
+      // CGAL vertex, then also insert the CGAL neighbours, otherwise
+      // just skip that step.
+      if (index != _supervertex[index].vertex->info().val()) continue;
+    } 
+
     // have a circulators that will go round the Voronoi neighbours of
     // _supervertex[index1].vertex
     Vertex_circulator vc = _TR.incident_vertices(_supervertex[index].vertex);
     Vertex_circulator done = vc;
-    do  {
-      // if a neighbouring vertex not the infinite vertex, then add it
-      // to our union of neighbouring vertices.
-      if (_verbose) cout << "examining " << vc->info().val() << endl;
-      if (vc->info().val() != INFINITE_VERTEX) {
-	// NB: from it=1 onwards occasionally it might already have
-	// been inserted -- but double insertion still leaves only one
-	// copy in the set, so there's no problem
-	NeighbourUnion.insert(vc->info().val());
-	if (_verbose) cout << "inserted " << vc->info().val() << endl;
-      } 
-    } while (++vc != done);
+    if (vc != NULL){ // a safety check in case there is no Voronoi
+		     // neighbour (which may happen e.g. if we just
+		     // have a bunch of coincident points)
+      do  {
+	// if a neighbouring vertex is not the infinite vertex, then add it
+	// to our union of neighbouring vertices.
+	if (_verbose) cout << "examining " << vc->info().val() << endl;
+	if (vc->info().val() != INFINITE_VERTEX) {
+	  // NB: from it=1 onwards occasionally it might already have
+	  // been inserted -- but double insertion still leaves only one
+	  // copy in the set, so there's no problem
+	  NeighbourUnion.insert(vc->info().val());
+	  if (_verbose) cout << "  inserted " << vc->info().val() << " to neighbours union" << endl;
+	} 
+      } while (++vc != done);
+    }
   }
   
@@ -160,9 +218,42 @@
   for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
     int index = indices_to_remove[ir];
+    if (_verbose) cout << "  removing " << index << endl;
 
     // NeighbourUnion should not contain the points to be removed
     // (because later we will assume they still exist).
     NeighbourUnion.erase(indices_to_remove[ir]);
-    
+
+    // first deal with  coincidences
+    if (_supervertex[index].coincidence != index){
+      int new_index = _supervertex[index].coincidence;
+
+      // if this is the point among the coiciding ones that "owns" the
+      // CGAL vertex we need to re-label the CGAL vertex so that it
+      // points to the coincident particle and set the current one to
+      // NULL
+      //
+      // This can be done only on the first point as they all share
+      // the same value
+      //
+      // Note that this has to be done before the following step since
+      // it will alter the coincidence information
+      if (index == _supervertex[index].vertex->info().val())
+	_supervertex[new_index].vertex->info() = new_index;
+
+      // we need to browse the coincidences until we end the loop, at
+      // which point we reset the coincidence of the point that has
+      // the current one as a coincidence
+      while (_supervertex[new_index].coincidence != index)
+	new_index = _supervertex[new_index].coincidence;
+      _supervertex[new_index].coincidence = _supervertex[index].coincidence;
+
+      // remove the coincidence on the point being removed and mark it
+      // as removed
+      _supervertex[index].coincidence = index;
+      _supervertex[index].vertex = NULL;
+
+      continue;
+    }
+
     // points to be removed should also be eliminated from the
     // triangulation and the supervertex structure should be updated
@@ -192,5 +283,7 @@
   // of the neighbour union happens to be on the wrong side.
   Face_handle face;
-  if (indices_to_remove.size() > 0) {
+  //if (indices_to_remove.size() > 0) { // GS: use NeighbourUnion instead
+                                        //     (safe also in case of coincidences)
+  if (NeighbourUnion.size() > 0) {
     // face can only be found if there were points to remove in first place
     face = _TR.incident_faces(
@@ -204,6 +297,8 @@
     int index = _supervertex.size()-1;
     indices_added.push_back(index);
-
-    if (indices_to_remove.size() > 0) {
+    if (_verbose) cout << "  adding " << index << endl;
+
+    //if (indices_to_remove.size() > 0) {
+    if (NeighbourUnion.size() > 0) {
       // be careful of using face (for location hinting) only when it exists
       _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first, 
@@ -213,8 +308,34 @@
 						    points_to_add[ia].second));
     }
-    // we are not up to dealing with coincident vertices, so make 
-    // sure the user knows!
-    _CrashIfVertexPresent(_supervertex[index].vertex, index);
-    _supervertex[index].vertex->info() = index;
+
+    // check if this leads to a coincidence
+    int coinciding_index = _CheckIfVertexPresent(_supervertex[index].vertex, index);
+    if (coinciding_index == index){
+      // we need to associate an index to each vertex -- thus when we get
+      // a vertex (e.g. as a nearest neighbour) from CGAL, we will be
+      // able to figure out which particle it corresponded to.
+      _supervertex[index].vertex->info() = _supervertex[index].coincidence = index;
+    } else {
+      if (_verbose) cout << "  coinciding with vertex " << coinciding_index << endl;
+      // the new vertex points to an already existing one and we
+      // record the coincidence
+      //
+      // we also update the NN of the coinciding particle (to avoid
+      // having to loop over the list of coinciding neighbours later)
+      // This is done first as it allows us to check if this is a new
+      // coincidence or a coincidence added to a particle that was
+      // previously "alone"
+      _supervertex[coinciding_index].NNindex = index;
+      _supervertex[coinciding_index].NNdistance = 0.0;
+      indices_of_updated_neighbours.push_back(coinciding_index);
+
+      // Note that we must not only set the coincidence of the
+      // currently-added particle, the one it coincides with also
+      // needs be updated (taking into account that it might already
+      // coincide with another one)
+      _supervertex[index].coincidence = _supervertex[coinciding_index].coincidence; // handles cases with previous coincidences
+      _supervertex[coinciding_index].coincidence = index;
+
+    }
     
     // first find nearest neighbour of "newpoint" (shorthand for
@@ -227,4 +348,8 @@
     indices_of_updated_neighbours.push_back(index);
     _SetAndUpdateNearest(index, indices_of_updated_neighbours);
+
+    //cout << "Added: " << setprecision(20) << " (" 
+    //     << points_to_add[ia].first << "," << points_to_add[ia].second
+    //     << ") with index " << index << endl;
   }
 
@@ -251,11 +376,77 @@
   }
 
+  if (_verbose) cout << "Leaving  DnnPlane::RemoveAndAddPoints" << endl;
 }
-
 
 //----------------------------------------------------------------------
 /// Determines the index and distance of the nearest neighbour to 
 /// point j and puts the information into the _supervertex entry for j.
-void DnnPlane::_SetNearest (const int & j) {
+void DnnPlane::_SetNearest (const int j) {
+  // first deal with the cases where we have a coincidence
+  if (_supervertex[j].coincidence != j){
+    _supervertex[j].NNindex = _supervertex[j].coincidence;
+    _supervertex[j].NNdistance = 0.0;
+    return;
+  }
+
+  // The code below entirely uses CGAL distance comparisons to compute
+  // the nearest neighbour. It has the mais drawback to induice a
+  // 10-20% time penalty so we switched to our own comparison (which
+  // only turns to CGAL for dangerous situations)
+  //
+  //  Vertex_handle current = _supervertex[j].vertex;
+  //  Vertex_circulator vc = _TR.incident_vertices(current);
+  //  Vertex_circulator done = vc;
+  //  Vertex_handle nearest = _TR.infinite_vertex();
+  //  double mindist = HUGE_DOUBLE;
+  //
+  //   // when there is only one finite point left in the triangulation, 
+  //   // there are no triangles. Presumably this is why voronoi returns
+  //   // NULL for the incident vertex circulator. Check if this is
+  //   // happening before circulating over it... (Otherwise it crashes
+  //   // when looking for neighbours of last point)
+  //   if (vc != NULL){
+  //     // initialise the nearest vertex handle to the first incident
+  //     // vertex that is not INFINITE_VERTEX
+  //     while (vc->info().val() == INFINITE_VERTEX){
+  //       vc++;
+  //       if (vc==done) break; // if vc==done, then INFINITE_VERTEX is the
+  // 			   // only element in the neighbourhood
+  //     }
+  // 
+  //     // if there is just the infinite vertex, we have vc->info().val()
+  //     // == INFINITE_VERTEX and nothing has to be done
+  //     // otherwise, use the current vc as an initialisation
+  //     if (vc->info().val() != INFINITE_VERTEX){
+  //       nearest = vc; // initialisation to the first non-infinite vertex
+  // 
+  //       // and loop over the following ones
+  //       while (++vc != done){
+  // 	// we should not compare with the infinite vertex
+  // 	if (vc->info().val() == INFINITE_VERTEX) continue;
+  // 
+  // 	if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl; 
+  // 	// use CGAL's distance comparison to check if 'vc' is closer to
+  // 	// 'current' than the nearest so far (we include the == case for
+  // 	// safety though it should not matter in this precise case)
+  // 	if (CGAL::compare_distance_to_point(current->point(), vc->point(), nearest->point())!=CGAL::LARGER){
+  // 	  nearest = vc;
+  // 	  if (_verbose) cout << "nearer";
+  // 	}
+  //       }
+  // 
+  //       // now compute the distance
+  //       //
+  //       // Note that since we're always using CGAL to compare distances
+  //       // (and never the distance computed using _euclid_distance) we
+  //       // should not worry about rounding errors in mindist
+  //       mindist = _euclid_distance(current->point(), nearest->point());
+  //     }
+  //   }
+  //
+  //  // set j's supervertex info about nearest neighbour
+  //  _supervertex[j].NNindex = nearest->info().val();
+  //  _supervertex[j].NNdistance = mindist;
+
   Vertex_handle current = _supervertex[j].vertex;
   Vertex_circulator vc = _TR.incident_vertices(current);
@@ -264,5 +455,5 @@
   double mindist = HUGE_DOUBLE; // change this to "HUGE" or max_double?
   Vertex_handle nearest = _TR.infinite_vertex();
-
+  
   // when there is only one finite point left in the triangulation, 
   // there are no triangles. Presumably this is why voronoi returns
@@ -274,14 +465,15 @@
       // find distance between j and its Voronoi neighbour (vc)
       if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
-      dist = _euclid_distance(current->point(), vc->point());
+
       // check if j is closer to vc than vc's currently registered
       // nearest neighbour (and update things if it is)
-      if (dist < mindist) {
-	mindist = dist; nearest = vc; 
-	if (_verbose) cout << "nearer ";
+      if (_is_closer_to(current->point(), vc->point(), nearest, dist, mindist)){
+	nearest = vc; 
+      	if (_verbose) cout << "nearer ";
       } 
       if (_verbose) cout << vc->point() << "; "<< dist << endl;
     }
   } while (++vc != done); // move on to next Voronoi neighbour
+  
   // set j's supervertex info about nearest neighbour
   _supervertex[j].NNindex = nearest->info().val();
@@ -291,5 +483,5 @@
 //----------------------------------------------------------------------
 /// Determines and stores the nearest neighbour of j, and where
-/// necessary update the nearest-neighbour info of Voronoi neighbours
+/// necessary updates the nearest-neighbour info of Voronoi neighbours
 /// of j;
 ///
@@ -304,10 +496,18 @@
 /// NB: note that we have _SetAndUpdateNearest as a completely
 ///     separate routine from _SetNearest because we want to
-///     use one single ciruclation over voronoi neighbours to find the
+///     use one single circulation over voronoi neighbours to find the
 ///     nearest neighbour and to update the voronoi neighbours if need
 ///     be.
 void DnnPlane::_SetAndUpdateNearest(
-			  const int & j, 
+			  const int j, 
 			  vector<int> & indices_of_updated_neighbours) {
+  //cout << "SetAndUpdateNearest for point " << j << endl;
+  // first deal with coincidences
+  if (_supervertex[j].coincidence != j){
+    _supervertex[j].NNindex = _supervertex[j].coincidence;
+    _supervertex[j].NNdistance = 0.0;
+    //cout << "  set to coinciding point " << _supervertex[j].coincidence << endl;
+    return;
+  }
 
   Vertex_handle current = _supervertex[j].vertex;
@@ -326,25 +526,45 @@
     if (vc->info().val() != INFINITE_VERTEX) {
       if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
-      // find distance between j and its Voronoi neighbour (vc)
-      dist = _euclid_distance(current->point(), vc->point());
+
       // update the mindist if we are closer than anything found so far
-      if (dist < mindist) {
-	mindist = dist; nearest = vc; 
-	if (_verbose) cout << "nearer ";
+      if (_is_closer_to(current->point(), vc->point(), nearest, dist, mindist)){
+	nearest = vc; 
+      	if (_verbose) cout << "nearer ";
       } 
+
       // find index corresponding to vc for easy manipulation
       int vcindx = vc->info().val();
       if (_verbose) cout << vc->point() << "; "<< dist << endl;
-      // check if j is closer to vc than vc's currently registered
-      // nearest neighbour (and update things if it is)
-      if (dist < _supervertex[vcindx].NNdistance) {
+
+      if (_is_closer_to_with_hint(vc->point(), current->point(), 
+				  _supervertex[_supervertex[vcindx].NNindex].vertex,
+				  dist, _supervertex[vcindx].NNdistance)){
 	if (_verbose) cout << vcindx << "'s NN becomes " << current->info().val() << endl;
-	_supervertex[vcindx].NNdistance = dist;
 	_supervertex[vcindx].NNindex = j;
 	indices_of_updated_neighbours.push_back(vcindx);
       }
+
+      // original code without the use of CGAL distance in potentially
+      // dangerous cases
+      //
+      // // check if j is closer to vc than vc's currently registered
+      // // nearest neighbour (and update things if it is)
+      // //
+      // // GS: originally, the distance test below was a strict <. It
+      // //     has to be <= because if the two distances are ==, it is
+      // //     possible that the old NN is no longer connected to vc in
+      // //     the triangulation, and we are sure that the newly
+      // //     inserted point (j) is (since we loop over j's
+      // //     neighbouring points in the triangulation).
+      // if (dist <= _supervertex[vcindx].NNdistance) {
+      // 	if (_verbose) cout << vcindx << "'s NN becomes " << current->info().val() << endl;
+      // 	_supervertex[vcindx].NNdistance = dist;
+      // 	_supervertex[vcindx].NNindex = j;
+      // 	indices_of_updated_neighbours.push_back(vcindx);
+      // }
     }
   } while (++vc != done); // move on to next Voronoi neighbour
   // set j's supervertex info about nearest neighbour
+  //cout << "  set to point " << nearest->info().val() << endl;
   _supervertex[j].NNindex = nearest->info().val();
   _supervertex[j].NNdistance = mindist;
