Context Navigation

-              r5b5a56b
+              r35cdc46
 //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
 //
 //----------------------------------------------------------------------
 …
 //
 //  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,
 …
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
 //ENDHEADER
+//FJENDHEADER
 …
 #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;
 …
        _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);
+  }
 …
 /// 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
 …
   // 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;
+}
 …
                           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
 …
   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);
+    }
+  }
 …
   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
 …
   // 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(
 …
     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,
 …
                                                     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
 …
     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;
+  }
 …
+  }
+  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);
 …
   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
 …
       // 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();
 …
 //----------------------------------------------------------------------
 /// 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;
 ///
 …
 /// 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;
 …
     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;

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 35cdc46 in git for external/fastjet/DnnPlane.cc

Legend:

external/fastjet/DnnPlane.cc

Download in other formats: