Changeset 23c605

Timestamp:

Aug 20, 2014, 1:06:16 PM (11 years ago)

Author:

Frederik Heber <heber@…>

Children:

1cde4e8

Parents:

1ae9aa

git-author:

Frederik Heber <heber@…> (06/30/14 09:35:42)

git-committer:

Frederik Heber <heber@…> (08/20/14 13:06:16)

Message:

SphericalPointDistribution is now working with bond degree weights.

• recurseMatching() now works on IndexTupleList_t.
• also rewrote calculatePairwiseDistances() and calculateErrorOfMatching().
• L1THRESHOLD in recurseMatching() moved over to class body.
• increased verbosity level of ...Matching() functions by one, added note on eventually chosen matching and why.
• we assert that bestL2 is not too large.
• FIX: calculateErrorOfMatching() did not use absolute value of gap for L1 error.
• TESTS: Regresssion test FragmentMolecule-cycles working again.

Files:

• src/Fragmentation/Exporters/SphericalPointDistribution.cpp (modified) (15 diffs)
• src/Fragmentation/Exporters/SphericalPointDistribution.hpp (modified) (4 diffs)
• src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp (modified) (1 diff)
• src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp (modified) (2 diffs)
• tests/regression/Fragmentation/FragmentMolecule/testsuite-fragmentation-fragment-molecule.at (modified) (1 diff)

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

@@ -49 +49 @@
 #include <limits>
 #include <list>
+#include <numeric>
 #include <vector>
 #include <map>
@@ -59 +60 @@
 // static entities
 const double SphericalPointDistribution::warn_amplitude = 1e-2;
+const double SphericalPointDistribution::L1THRESHOLD = 1e-2;
+const double SphericalPointDistribution::L2THRESHOLD = 2e-1;
 
 typedef std::vector<double> DistanceArray_t;
@@ -64 +67 @@
 // class generator: taken from www.cplusplus.com example std::generate
 struct c_unique {
-  int current;
+  unsigned int current;
   c_unique() {current=0;}
-  int operator()() {return current++;}
+  unsigned int operator()() {return current++;}
 } UniqueNumber;
 
-inline
-DistanceArray_t calculatePairwiseDistances(
-    const std::vector<Vector> &_points,
-    const SphericalPointDistribution::IndexList_t &_indices
-    )
-{
-  DistanceArray_t result;
-  for (SphericalPointDistribution::IndexList_t::const_iterator firstiter = _indices.begin();
-      firstiter != _indices.end(); ++firstiter) {
-    for (SphericalPointDistribution::IndexList_t::const_iterator seconditer = firstiter;
-        seconditer != _indices.end(); ++seconditer) {
-      if (firstiter == seconditer)
-        continue;
-      const double distance = (_points[*firstiter] - _points[*seconditer]).NormSquared();
-      result.push_back(distance);
-    }
-  }
-  return result;
-}
+struct c_unique_list {
+  unsigned int current;
+  c_unique_list() {current=0;}
+  std::list<unsigned int> operator()() {return std::list<unsigned int>(1, current++);}
+} UniqueNumberList;
 
 /** Calculate the center of a given set of points in \a _positions but only
@@ -107 +96 @@
 
   return Center;
+}
+
+inline
+DistanceArray_t calculatePairwiseDistances(
+    const SphericalPointDistribution::VectorArray_t &_points,
+    const SphericalPointDistribution::IndexTupleList_t &_indices
+    )
+{
+  DistanceArray_t result;
+  for (SphericalPointDistribution::IndexTupleList_t::const_iterator firstiter = _indices.begin();
+      firstiter != _indices.end(); ++firstiter) {
+
+    // calculate first center from possible tuple of indices
+    Vector FirstCenter;
+    ASSERT(!firstiter->empty(), "calculatePairwiseDistances() - there is an empty tuple.");
+    if (firstiter->size() == 1) {
+      FirstCenter = _points[*firstiter->begin()];
+    } else {
+      FirstCenter = calculateCenter( _points, *firstiter);
+      if (!FirstCenter.IsZero())
+        FirstCenter.Normalize();
+    }
+
+    for (SphericalPointDistribution::IndexTupleList_t::const_iterator seconditer = firstiter;
+        seconditer != _indices.end(); ++seconditer) {
+      if (firstiter == seconditer)
+        continue;
+
+      // calculate second center from possible tuple of indices
+      Vector SecondCenter;
+      ASSERT(!seconditer->empty(), "calculatePairwiseDistances() - there is an empty tuple.");
+      if (seconditer->size() == 1) {
+        SecondCenter = _points[*seconditer->begin()];
+      } else {
+        SecondCenter = calculateCenter( _points, *seconditer);
+        if (!SecondCenter.IsZero())
+          SecondCenter.Normalize();
+      }
+
+      // calculate distance between both centers
+      double distance = 2.; // greatest distance on surface of sphere with radius 1.
+      if ((!FirstCenter.IsZero()) && (!SecondCenter.IsZero()))
+        distance = (FirstCenter - SecondCenter).NormSquared();
+      result.push_back(distance);
+    }
+  }
+  return result;
 }
 
@@ -123 +159 @@
   Vector dreiBein(_oldPosition);
   dreiBein.VectorProduct(_RotationAxis);
+  ASSERT( !dreiBein.IsZero(), "determineSignOfRotation() - dreiBein is zero.");
   dreiBein.Normalize();
   const double sign =
       (dreiBein.ScalarProduct(_newPosition) < 0.) ? -1. : +1.;
   LOG(6, "DEBUG: oldCenter on plane is " << _oldPosition
-      << ", newCenter in plane is " << _newPosition
+      << ", newCenter on plane is " << _newPosition
       << ", and dreiBein is " << dreiBein);
   return sign;
@@ -166 +203 @@
  */
 std::pair<double, double> SphericalPointDistribution::calculateErrorOfMatching(
-    const std::vector<Vector> &_old,
-    const std::vector<Vector> &_new,
-    const IndexList_t &_Matching)
+    const VectorArray_t &_old,
+    const VectorArray_t &_new,
+    const IndexTupleList_t &_Matching)
 {
   std::pair<double, double> errors( std::make_pair( 0., 0. ) );
 
   if (_Matching.size() > 1) {
-    LOG(3, "INFO: Matching is " << _Matching);
+    LOG(5, "INFO: Matching is " << _Matching);
 
     // calculate all pair-wise distances
-    IndexList_t keys(_Matching.size());
-    std::generate (keys.begin(), keys.end(), UniqueNumber);
+    IndexTupleList_t keys(_old.size(), IndexList_t() );
+    std::generate (keys.begin(), keys.end(), UniqueNumberList);
+
     const DistanceArray_t firstdistances = calculatePairwiseDistances(_old, keys);
     const DistanceArray_t seconddistances = calculatePairwiseDistances(_new, _Matching);
@@ -187 +225 @@
     for (;(firstiter != firstdistances.end()) && (seconditer != seconddistances.end());
         ++firstiter, ++seconditer) {
-      const double gap = *firstiter - *seconditer;
+      const double gap = fabs(*firstiter - *seconditer);
       // L1 error
       if (errors.first < gap)
@@ -194 +232 @@
       errors.second += gap*gap;
     }
-  } else
-    ELOG(3, "calculateErrorOfMatching() - Given matching's size is less than 2.");
-  LOG(3, "INFO: Resulting errors for matching (L1, L2): "
+  } else {
+    // check whether we have any zero centers: Combining points on new sphere may result
+    // in zero centers
+    for (SphericalPointDistribution::IndexTupleList_t::const_iterator iter = _Matching.begin();
+        iter != _Matching.end(); ++iter) {
+      if ((iter->size() != 1) && (calculateCenter( _new, *iter).IsZero())) {
+        errors.first = 2.;
+        errors.second = 2.;
+      }
+    }
+  }
+  LOG(4, "INFO: Resulting errors for matching (L1, L2): "
       << errors.first << "," << errors.second << ".");
 
@@ -232 +279 @@
  * \param _matching current matching being build up
  * \param _indices contains still available indices
+ * \param _remainingweights current weights to fill (each weight a place)
+ * \param _remainiter iterator over the weights, indicating the current position we match
  * \param _matchingsize
  */
 void SphericalPointDistribution::recurseMatchings(
     MatchingControlStructure &_MCS,
-    IndexList_t &_matching,
+    IndexTupleList_t &_matching,
     IndexList_t _indices,
-    unsigned int _matchingsize)
-{
-  LOG(4, "DEBUG: Recursing with current matching " << _matching
+    WeightsArray_t &_remainingweights,
+    WeightsArray_t::iterator _remainiter,
+    const unsigned int _matchingsize
+    )
+{
+  LOG(5, "DEBUG: Recursing with current matching " << _matching
       << ", remaining indices " << _indices
-      << ", and sought size " << _matching.size()+_matchingsize);
-  //!> threshold for L1 error below which matching is immediately acceptable
-  const double L1THRESHOLD = 1e-2;
+      << ", and remaining weights " << _matchingsize);
   if (!_MCS.foundflag) {
-    LOG(4, "DEBUG: Current matching has size " << _matching.size() << ", places left " << _matchingsize);
+    LOG(5, "DEBUG: Current matching has size " << _matching.size() << ", places left " << _matchingsize);
     if (_matchingsize > 0) {
       // go through all indices
       for (IndexList_t::iterator iter = _indices.begin();
           (iter != _indices.end()) && (!_MCS.foundflag);) {
+        // check whether we can stay in the current bin or have to move on to next one
+        if (*_remainiter == 0) {
+          // we need to move on
+          if (_remainiter != _remainingweights.end()) {
+            ++_remainiter;
+          } else {
+            // as we check _matchingsize > 0 this should be impossible
+            ASSERT( 0, "recurseMatchings() - we must not come to this position.");
+          }
+        }
+        // advance in matching to same position
+        const size_t OldIndex = std::distance(_remainingweights.begin(), _remainiter);
+        while (_matching.size() <= OldIndex) { // add empty lists of new bin is opened
+          LOG(6, "DEBUG: Extending _matching.");
+          _matching.push_back( IndexList_t() );
+        }
+        IndexTupleList_t::iterator filliniter = _matching.begin();
+        std::advance(filliniter, OldIndex);
         // add index to matching
-        _matching.push_back(*iter);
-        LOG(5, "DEBUG: Adding " << *iter << " to matching.");
+        filliniter->push_back(*iter);
+        --(*_remainiter);
+        LOG(6, "DEBUG: Adding " << *iter << " to matching at " << OldIndex << ".");
         // remove index but keep iterator to position (is the next to erase element)
         IndexList_t::iterator backupiter = _indices.erase(iter);
         // recurse with decreased _matchingsize
-        recurseMatchings(_MCS, _matching, _indices, _matchingsize-1);
+        recurseMatchings(_MCS, _matching, _indices, _remainingweights, _remainiter, _matchingsize-1);
         // re-add chosen index and reset index to new position
-        _indices.insert(backupiter, _matching.back());
+        _indices.insert(backupiter, filliniter->back());
         iter = backupiter;
         // remove index from _matching to make space for the next one
-        _matching.pop_back();
+        filliniter->pop_back();
+        ++(*_remainiter);
       }
       // gone through all indices then exit recursion
@@ -268 +338 @@
         _MCS.foundflag = true;
     } else {
-      LOG(3, "INFO: Found matching " << _matching);
+      LOG(4, "INFO: Found matching " << _matching);
       // calculate errors
       std::pair<double, double> errors = calculateErrorOfMatching(
@@ -309 +379 @@
   MCS.foundflag = false;
   MCS.bestL2 = std::numeric_limits<double>::max();
+  // transform lists into arrays
   for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
-      iter != _polygon.end(); ++iter)
+      iter != _polygon.end(); ++iter) {
     MCS.oldpoints.push_back(iter->first);
+    MCS.weights.push_back(iter->second);
+  }
   MCS.newpoints.insert(MCS.newpoints.begin(), _newpolygon.begin(),_newpolygon.end() );
 
@@ -319 +392 @@
     IndexList_t indices(_newpolygon.size());
     std::generate(indices.begin(), indices.end(), UniqueNumber);
-    IndexList_t matching;
+    IndexTupleList_t matching;
 
     // walk through all matchings
-    const unsigned int matchingsize = _polygon.size();
-    ASSERT( matchingsize <= indices.size(),
-        "SphericalPointDistribution::matchSphericalPointDistributions() - not enough new points to choose for matching to old ones.");
-    recurseMatchings(MCS, matching, indices, matchingsize);
+    WeightsArray_t remainingweights = MCS.weights;
+    unsigned int placesleft = std::accumulate(remainingweights.begin(), remainingweights.end(), 0);
+    recurseMatchings(MCS, matching, indices, remainingweights, remainingweights.begin(), placesleft);
+  }
+  if (MCS.foundflag)
+    LOG(3, "Found a best matching beneath L1 threshold of " << L1THRESHOLD);
+  else {
+    if (MCS.bestL2 < warn_amplitude)
+      LOG(3, "Picking matching is " << MCS.bestmatching << " with best L2 error of "
+          << MCS.bestL2);
+    else if (MCS.bestL2 < L2THRESHOLD)
+      ELOG(2, "Picking matching is " << MCS.bestmatching
+          << " with rather large L2 error of " << MCS.bestL2);
+    else
+      ASSERT(0, "findBestMatching() - matching "+toString(MCS.bestmatching)
+          +" has L2 error of "+toString(MCS.bestL2)+" that is too large.");
   }
 
   // combine multiple points and create simple IndexList from IndexTupleList
-  IndexTupleList_t bestmatching;
-  for (IndexList_t::const_iterator iter = MCS.bestmatching.begin();
-      iter != MCS.bestmatching.end(); ++iter)
-    bestmatching.push_back(IndexList_t(1, *iter));
   const SphericalPointDistribution::IndexList_t IndexList =
-      joinPoints(_newpolygon, MCS.newpoints, bestmatching);
+      joinPoints(_newpolygon, MCS.newpoints, MCS.bestmatching);
 
   return IndexList;
@@ -363 +444 @@
       Center.Normalize();
       _newpolygon.push_back(Center);
-      LOG(5, "DEBUG: Combining " << tupleiter->size() << "points to weighted center "
+      LOG(5, "DEBUG: Combining " << tupleiter->size() << " points to weighted center "
           << Center << " with new index " << UniqueIndex);
       // mark for removal
@@ -600 +681 @@
       iter != _polygon.end(); ++iter)
     remainingold.push_back(iter->first);
-  VectorArray_t remainingnew(_newpolygon.begin(), _newpolygon.end());
   LOG(2, "INFO: Matching old polygon " << _polygon
       << " with new polygon " << _newpolygon);
@@ -613 +693 @@
     IndexList_t bestmatching = findBestMatching(_polygon, _newpolygon);
     LOG(2, "INFO: Best matching is " << bestmatching);
+
+    // _newpolygon has changed, so now convert to array
+    VectorArray_t remainingnew(_newpolygon.begin(), _newpolygon.end());
 
     // determine rotation angles to align the two point distributions with

src/Fragmentation/Exporters/SphericalPointDistribution.hpp

@@ -89 +89 @@
   //!> precalculated value for root of 3
   const double SQRT_3;
+  //!> threshold for L1 error below which matching is immediately acceptable
+  static const double L1THRESHOLD;
+  //!> threshold for L2 error below which matching is acceptable
+  static const double L2THRESHOLD;
 
   //!> typedef for a full rotation specification consisting of axis and angle.
@@ -103 +107 @@
   //!> typedef for a Vector of positions with weights
   typedef std::vector< std::pair<Vector, int> > WeightedVectorArray_t;
+  //!> typedef for a vector of degrees (or integral weights)
+  typedef std::vector<unsigned int> WeightsArray_t;
 
   //!> amplitude up to which deviations in checks of rotations are tolerated
@@ -112 +118 @@
 
   static std::pair<double, double> calculateErrorOfMatching(
-      const std::vector<Vector> &_old,
-      const std::vector<Vector> &_new,
-      const IndexList_t &_Matching);
+      const VectorArray_t &_old,
+      const VectorArray_t &_new,
+      const IndexTupleList_t &_Matching);
 
   static Polygon_t removeMatchingPoints(
@@ -124 +130 @@
     bool foundflag;
     double bestL2;
-    IndexList_t bestmatching;
+    IndexTupleList_t bestmatching;
     VectorArray_t oldpoints;
     VectorArray_t newpoints;
+    WeightsArray_t weights;
   };
 
   static void recurseMatchings(
       MatchingControlStructure &_MCS,
-      IndexList_t &_matching,
+      IndexTupleList_t &_matching,
       IndexList_t _indices,
-      unsigned int _matchingsize);
+      WeightsArray_t &_remainingweights,
+      WeightsArray_t::iterator _remainiter,
+      const unsigned int _matchingsize
+      );
 
   static IndexList_t findBestMatching(

src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp

@@ -638 +638 @@
 }
 
+/** UnitTest for matchSphericalPointDistributions() with four points and weights
+ * not all equal to one.
+ */
+void SphericalPointDistributionTest::matchSphericalPointDistributionsTest_multiple()
+{
+  SphericalPointDistribution SPD(1.);
+
+  // test with four points: one point having weight of two
+  {
+    SphericalPointDistribution::WeightedPolygon_t polygon;
+    polygon += std::make_pair( Vector(1.,0.,0.), 2);
+    SphericalPointDistribution::Polygon_t newpolygon =
+        SPD.get<4>();
+    SphericalPointDistribution::Polygon_t expected;
+    expected += Vector(-0.5773502691896,-5.551115123126e-17,0.8164965809277);
+    expected += Vector(-0.5773502691896,-5.551115123126e-17,-0.8164965809277);
+    SphericalPointDistribution::Polygon_t remaining =
+        SphericalPointDistribution::matchSphericalPointDistributions(
+            polygon,
+            newpolygon);
+//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
+    CPPUNIT_ASSERT_EQUAL( expected, remaining );
+  }
+
+  // test with five points: one point having weight of two
+  {
+    SphericalPointDistribution::WeightedPolygon_t polygon;
+    polygon += std::make_pair( Vector(1.,0.,0.), 2);
+    SphericalPointDistribution::Polygon_t newpolygon =
+        SPD.get<5>();
+    SphericalPointDistribution::Polygon_t expected;
+    expected += Vector(-0.7071067811865,0.7071067811865,0);
+    expected += Vector(-0.3535533905933,-0.3535533905933,0.8660254037844);
+    expected += Vector(-0.3535533905933,-0.3535533905933,-0.8660254037844);
+    SphericalPointDistribution::Polygon_t remaining =
+        SphericalPointDistribution::matchSphericalPointDistributions(
+            polygon,
+            newpolygon);
+//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
+    CPPUNIT_ASSERT_EQUAL( expected, remaining );
+  }
+
+
+  // test with five points: one point having weight of two, one weight of one
+  {
+    SphericalPointDistribution::WeightedPolygon_t polygon;
+    polygon += std::make_pair( Vector(M_SQRT1_2,M_SQRT1_2,0.), 2);
+    polygon += std::make_pair( Vector(-1.,0.,0.), 1);
+    SphericalPointDistribution::Polygon_t newpolygon =
+        SPD.get<5>();
+    SphericalPointDistribution::Polygon_t expected;
+    expected += Vector(0.3535533786708,-0.3535533955317,-0.8660254066357);
+    expected += Vector(0.3535534025157,-0.3535533856548,0.8660254009332);
+    SphericalPointDistribution::Polygon_t remaining =
+        SphericalPointDistribution::matchSphericalPointDistributions(
+            polygon,
+            newpolygon);
+//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
+    CPPUNIT_ASSERT_EQUAL( expected, remaining );
+  }
+
+  // test with six points: two points each having weight of two
+  {
+    SphericalPointDistribution::WeightedPolygon_t polygon;
+    polygon += std::make_pair( Vector(M_SQRT1_2,-M_SQRT1_2,0.), 2);
+    polygon += std::make_pair( Vector(-M_SQRT1_2,M_SQRT1_2,0.), 2);
+    SphericalPointDistribution::Polygon_t newpolygon =
+        SPD.get<6>();
+    SphericalPointDistribution::Polygon_t expected;
+    expected += Vector(0.,0.,-1.);
+    expected += Vector(0.,0.,1.);
+    SphericalPointDistribution::Polygon_t remaining =
+        SphericalPointDistribution::matchSphericalPointDistributions(
+            polygon,
+            newpolygon);
+//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
+    CPPUNIT_ASSERT_EQUAL( expected, remaining );
+  }
+}
+
 /** UnitTest for matchSphericalPointDistributions() with five points
  */

src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp

@@ -31 +31 @@
     CPPUNIT_TEST ( matchSphericalPointDistributionsTest_7 );
     CPPUNIT_TEST ( matchSphericalPointDistributionsTest_8 );
+    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_multiple );
     CPPUNIT_TEST_SUITE_END();
 
@@ -45 +46 @@
       void matchSphericalPointDistributionsTest_7();
       void matchSphericalPointDistributionsTest_8();
+      void matchSphericalPointDistributionsTest_multiple();
 
 private:

tests/regression/Fragmentation/FragmentMolecule/testsuite-fragmentation-fragment-molecule.at

@@ -37 +37 @@
 AT_SETUP([Fragmentation - Fragmentation with cycles])
 AT_KEYWORDS([fragmentation fragment-molecule cycle])
-AT_XFAIL_IF([/bin/true])
 
 file=benzene.pdb