Changeset 450adf for src/Fragmentation

Timestamp:

Sep 12, 2016, 11:48:36 PM (8 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_StructOpt_integration_tests, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, GeometryObjects, Gui_displays_atomic_force_velocity, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, Ubuntu_1604_changes, stable

Children:

653cea

Parents:

2199c2

git-author:

Frederik Heber <heber@…> (06/29/14 21:20:49)

git-committer:

Frederik Heber <heber@…> (09/12/16 23:48:36)

Message:

Extracted joinPoints() function to make it accessible to tests.

• added unit test.
• moved some function definitions around and added documentation.

Location:

src/Fragmentation/Exporters

Files:

• SphericalPointDistribution.cpp (modified) (7 diffs)
• SphericalPointDistribution.hpp (modified) (5 diffs)
• unittests/SphericalPointDistributionUnitTest.cpp (modified) (1 diff)
• unittests/SphericalPointDistributionUnitTest.hpp (modified) (2 diffs)

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

@@ -63 +63 @@
 typedef std::vector<double> DistanceArray_t;
 
+// class generator: taken from www.cplusplus.com example std::generate
+struct c_unique {
+  int current;
+  c_unique() {current=0;}
+  int operator()() {return current++;}
+} UniqueNumber;
+
 inline
 DistanceArray_t calculatePairwiseDistances(
@@ -83 +90 @@
 }
 
-// class generator: taken from www.cplusplus.com example std::generate
-struct c_unique {
-  int current;
-  c_unique() {current=0;}
-  int operator()() {return current++;}
-} UniqueNumber;
-
+/** Calculate the center of a given set of points in \a _positions but only
+ * for those indicated by \a _indices.
+ *
+ */
+inline
+Vector calculateCenter(
+    const SphericalPointDistribution::VectorArray_t &_positions,
+    const SphericalPointDistribution::IndexList_t &_indices)
+{
+  Vector Center;
+  Center.Zero();
+  for (SphericalPointDistribution::IndexList_t::const_iterator iter = _indices.begin();
+      iter != _indices.end(); ++iter)
+    Center += _positions[*iter];
+  if (!_indices.empty())
+    Center *= 1./(double)_indices.size();
+
+  return Center;
+}
+
+/** Decides by an orthonormal third vector whether the sign of the rotation
+ * angle should be negative or positive.
+ *
+ * \return -1 or 1
+ */
+inline
+double determineSignOfRotation(
+    const Vector &_oldPosition,
+    const Vector &_newPosition,
+    const Vector &_RotationAxis
+    )
+{
+  Vector dreiBein(_oldPosition);
+  dreiBein.VectorProduct(_RotationAxis);
+  dreiBein.Normalize();
+  const double sign =
+      (dreiBein.ScalarProduct(_newPosition) < 0.) ? -1. : +1.;
+  LOG(6, "DEBUG: oldCenter on plane is " << _oldPosition
+      << ", newCenter in plane is " << _newPosition
+      << ", and dreiBein is " << dreiBein);
+  return sign;
+}
+
+/** Convenience function to recalculate old and new center for determining plane
+ * rotation.
+ */
+inline
+void calculateOldAndNewCenters(
+    Vector &_oldCenter,
+    Vector &_newCenter,
+    const SphericalPointDistribution::VectorArray_t &_referencepositions,
+    const SphericalPointDistribution::VectorArray_t &_currentpositions,
+    const SphericalPointDistribution::IndexList_t &_bestmatching)
+{
+  const size_t NumberIds = std::min(_bestmatching.size(), (size_t)3);
+  SphericalPointDistribution::IndexList_t continuousIds(NumberIds, -1);
+  std::generate(continuousIds.begin(), continuousIds.end(), UniqueNumber);
+  _oldCenter = calculateCenter(_referencepositions, continuousIds);
+  // C++11 defines a copy_n function ...
+  SphericalPointDistribution::IndexList_t::const_iterator enditer = _bestmatching.begin();
+  std::advance(enditer, NumberIds);
+  SphericalPointDistribution::IndexList_t firstbestmatchingIds(NumberIds, -1);
+  std::copy(_bestmatching.begin(), enditer, firstbestmatchingIds.begin());
+  _newCenter = calculateCenter( _currentpositions, firstbestmatchingIds);
+}
 /** Returns squared L2 error of the given \a _Matching.
  *
@@ -219 +284 @@
 }
 
-/** Decides by an orthonormal third vector whether the sign of the rotation
- * angle should be negative or positive.
- *
- * \return -1 or 1
- */
-inline
-double determineSignOfRotation(
-    const Vector &_oldPosition,
-    const Vector &_newPosition,
-    const Vector &_RotationAxis
-    )
-{
-  Vector dreiBein(_oldPosition);
-  dreiBein.VectorProduct(_RotationAxis);
-  dreiBein.Normalize();
-  const double sign =
-      (dreiBein.ScalarProduct(_newPosition) < 0.) ? -1. : +1.;
-  LOG(6, "DEBUG: oldCenter on plane is " << _oldPosition
-      << ", newCenter in plane is " << _newPosition
-      << ", and dreiBein is " << dreiBein);
-  return sign;
-}
-
 /** Finds combinatorially the best matching between points in \a _polygon
  * and \a _newpolygon.
@@ -247 +289 @@
  * We find the matching with the smallest L2 error, where we break when we stumble
  * upon a matching with zero error.
+ *
+ * As points in \a _polygon may be have a weight greater 1 we have to match it to
+ * multiple points in \a _newpolygon. Eventually, these multiple points are combined
+ * for their center of weight, which is the only thing follow-up function see of
+ * these multiple points. Hence, we actually modify \a _newpolygon in the process
+ * such that the returned IndexList_t indicates a bijective mapping in the end.
  *
  * \sa recurseMatchings() for going through all matchings
@@ -255 +303 @@
  */
 SphericalPointDistribution::IndexList_t SphericalPointDistribution::findBestMatching(
-    const SphericalPointDistribution::WeightedPolygon_t &_polygon,
-    const SphericalPointDistribution::Polygon_t &_newpolygon
+    const WeightedPolygon_t &_polygon,
+    Polygon_t &_newpolygon
     )
 {
@@ -280 +328 @@
     recurseMatchings(MCS, matching, indices, matchingsize);
   }
-  return MCS.bestmatching;
-}
-
-inline
-Vector calculateCenter(
-    const SphericalPointDistribution::VectorArray_t &_positions,
-    const SphericalPointDistribution::IndexList_t &_indices)
-{
-  Vector Center;
-  Center.Zero();
-  for (SphericalPointDistribution::IndexList_t::const_iterator iter = _indices.begin();
-      iter != _indices.end(); ++iter)
-    Center += _positions[*iter];
-  if (!_indices.empty())
-    Center *= 1./(double)_indices.size();
-
-  return Center;
-}
-
-inline
-void calculateOldAndNewCenters(
-    Vector &_oldCenter,
-    Vector &_newCenter,
-    const SphericalPointDistribution::VectorArray_t &_referencepositions,
-    const SphericalPointDistribution::VectorArray_t &_currentpositions,
-    const SphericalPointDistribution::IndexList_t &_bestmatching)
-{
-  const size_t NumberIds = std::min(_bestmatching.size(), (size_t)3);
-  SphericalPointDistribution::IndexList_t continuousIds(NumberIds, -1);
-  std::generate(continuousIds.begin(), continuousIds.end(), UniqueNumber);
-  _oldCenter = calculateCenter(_referencepositions, continuousIds);
-  // C++11 defines a copy_n function ...
-  SphericalPointDistribution::IndexList_t::const_iterator enditer = _bestmatching.begin();
-  std::advance(enditer, NumberIds);
-  SphericalPointDistribution::IndexList_t firstbestmatchingIds(NumberIds, -1);
-  std::copy(_bestmatching.begin(), enditer, firstbestmatchingIds.begin());
-  _newCenter = calculateCenter( _currentpositions, firstbestmatchingIds);
+
+  // 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);
+
+  return IndexList;
+}
+
+SphericalPointDistribution::IndexList_t SphericalPointDistribution::joinPoints(
+    Polygon_t &_newpolygon,
+    const VectorArray_t &_newpoints,
+    const IndexTupleList_t &_bestmatching
+    )
+{
+  // combine all multiple points
+  IndexList_t IndexList;
+  IndexArray_t removalpoints;
+  unsigned int UniqueIndex = _newpolygon.size(); // all indices up to size are used right now
+  VectorArray_t newCenters;
+  newCenters.reserve(_bestmatching.size());
+  for (IndexTupleList_t::const_iterator tupleiter = _bestmatching.begin();
+      tupleiter != _bestmatching.end(); ++tupleiter) {
+    ASSERT (tupleiter->size() > 0,
+        "findBestMatching() - encountered tuple in bestmatching with size 0.");
+    if (tupleiter->size() == 1) {
+      // add point and index
+      IndexList.push_back(*tupleiter->begin());
+    } else {
+      // combine into weighted and normalized center
+      Vector Center = calculateCenter(_newpoints, *tupleiter);
+      Center.Normalize();
+      _newpolygon.push_back(Center);
+      LOG(5, "DEBUG: Combining " << tupleiter->size() << "points to weighted center "
+          << Center << " with new index " << UniqueIndex);
+      // mark for removal
+      removalpoints.insert(removalpoints.end(), tupleiter->begin(), tupleiter->end());
+      // add new index
+      IndexList.push_back(UniqueIndex++);
+    }
+  }
+  // IndexList is now our new bestmatching (that is bijective)
+  LOG(4, "DEBUG: Our new bijective IndexList reads as " << IndexList);
+
+  // modifying _newpolygon: remove all points in removalpoints, add those in newCenters
+  Polygon_t allnewpoints = _newpolygon;
+  {
+    _newpolygon.clear();
+    std::sort(removalpoints.begin(), removalpoints.end());
+    size_t i = 0;
+    IndexArray_t::const_iterator removeiter = removalpoints.begin();
+    for (Polygon_t::iterator iter = allnewpoints.begin();
+        iter != allnewpoints.end(); ++iter, ++i) {
+      if ((removeiter != removalpoints.end()) && (i == *removeiter)) {
+        // don't add, go to next remove index
+        ++removeiter;
+      } else {
+        // otherwise add points
+        _newpolygon.push_back(*iter);
+      }
+    }
+  }
+  LOG(4, "DEBUG: The polygon with recentered points removed is " << _newpolygon);
+
+  // map IndexList to new shrinked _newpolygon
+  typedef std::set<unsigned int> IndexSet_t;
+  IndexSet_t SortedIndexList(IndexList.begin(), IndexList.end());
+  IndexList.clear();
+  {
+    size_t offset = 0;
+    IndexSet_t::const_iterator listiter = SortedIndexList.begin();
+    IndexArray_t::const_iterator removeiter = removalpoints.begin();
+    for (size_t i = 0; i < allnewpoints.size(); ++i) {
+      if ((removeiter != removalpoints.end()) && (i == *removeiter)) {
+        ++offset;
+        ++removeiter;
+      } else if ((listiter != SortedIndexList.end()) && (i == *listiter)) {
+        IndexList.push_back(*listiter - offset);
+        ++listiter;
+      }
+    }
+  }
+  LOG(4, "DEBUG: Our new bijective IndexList corrected for removed points reads as "
+      << IndexList);
+
+  return IndexList;
 }
 
@@ -497 +597 @@
 SphericalPointDistribution::matchSphericalPointDistributions(
     const SphericalPointDistribution::WeightedPolygon_t &_polygon,
-    const SphericalPointDistribution::Polygon_t &_newpolygon
+    SphericalPointDistribution::Polygon_t &_newpolygon
     )
 {

src/Fragmentation/Exporters/SphericalPointDistribution.hpp

@@ -21 +21 @@
 
 #include "LinearAlgebra/Vector.hpp"
+
+class SphericalPointDistributionTest;
 
 /** contains getters for the VSEPR model for specific number of electrons.
@@ -72 +74 @@
   static Polygon_t matchSphericalPointDistributions(
       const WeightedPolygon_t &_polygon,
-      const Polygon_t &_newpolygon
+      Polygon_t &_newpolygon
       );
 
@@ -80 +82 @@
   static const double SQRT_3;
 
+  //!> typedef for a full rotation specification consisting of axis and angle.
   typedef std::pair<Vector, double> Rotation_t;
 
+  //!> typedef for a list of indices (of points in a polygon)
   typedef std::list<unsigned int> IndexList_t;
+  //!> typedef enumerating possibly multiple points accumulated as one point
+  typedef std::list< IndexList_t > IndexTupleList_t;
+  //!> typedef for a vector of indices
   typedef std::vector<unsigned int> IndexArray_t;
+  //!> typedef for a Vector of positions
   typedef std::vector<Vector> VectorArray_t;
+  //!> typedef for a Vector of positions with weights
+  typedef std::vector< std::pair<Vector, int> > WeightedVectorArray_t;
 
   //!> amplitude up to which deviations in checks of rotations are tolerated
@@ -90 +100 @@
 
 private:
+  //!> grant unit tests access to private parts
+  friend class SphericalPointDistributionTest;
+
   static std::pair<double, double> calculateErrorOfMatching(
       const std::vector<Vector> &_old,
@@ -116 +129 @@
   static IndexList_t findBestMatching(
       const WeightedPolygon_t &_polygon,
-      const Polygon_t &_newpolygon
+      Polygon_t &_newpolygon
+      );
+
+  static IndexList_t joinPoints(
+      Polygon_t &_newpolygon,
+      const VectorArray_t &_newpoints,
+      const IndexTupleList_t &_bestmatching
       );
 

src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp

@@ -245 +245 @@
 }
 
+/** UnitTest for joinPoints()
+ */
+void SphericalPointDistributionTest::joinPointsTest()
+{
+  // test with simple configuration of three points
+  {
+    SphericalPointDistribution::Polygon_t newpolygon;
+    newpolygon += Vector(1.,0.,0.);
+    newpolygon += Vector(0.,1.,0.);
+    newpolygon += Vector(0.,0.,1.);
+    SphericalPointDistribution::Polygon_t expectedpolygon = newpolygon;
+    SphericalPointDistribution::IndexTupleList_t matching;
+    matching += SphericalPointDistribution::IndexList_t(1,0);
+    matching += SphericalPointDistribution::IndexList_t(1,1);
+    matching += SphericalPointDistribution::IndexList_t(1,2);
+    SphericalPointDistribution::IndexList_t IndexList =
+        SphericalPointDistribution::joinPoints(
+            newpolygon,
+            SphericalPointDistribution::VectorArray_t(newpolygon.begin(), newpolygon.end()),
+            matching);
+    SphericalPointDistribution::IndexList_t expected;
+    expected += 0,1,2;
+    CPPUNIT_ASSERT_EQUAL( expected, IndexList );
+    CPPUNIT_ASSERT_EQUAL( expectedpolygon, newpolygon );
+  }
+
+  // test with simple configuration of three points, only two are picked
+  {
+    SphericalPointDistribution::Polygon_t newpolygon;
+    newpolygon += Vector(1.,0.,0.);
+    newpolygon += Vector(0.,1.,0.);
+    newpolygon += Vector(0.,0.,1.);
+    SphericalPointDistribution::Polygon_t expectedpolygon = newpolygon;
+    SphericalPointDistribution::IndexTupleList_t matching;
+    matching += SphericalPointDistribution::IndexList_t(1,1);
+    matching += SphericalPointDistribution::IndexList_t(1,2);
+    SphericalPointDistribution::IndexList_t IndexList =
+        SphericalPointDistribution::joinPoints(
+            newpolygon,
+            SphericalPointDistribution::VectorArray_t(newpolygon.begin(), newpolygon.end()),
+            matching);
+    SphericalPointDistribution::IndexList_t expected;
+    expected += 1,2;
+    CPPUNIT_ASSERT_EQUAL( expected, IndexList );
+    CPPUNIT_ASSERT_EQUAL( expectedpolygon, newpolygon );
+  }
+
+  // test with simple configuration of three points, two are joined
+  {
+    SphericalPointDistribution::Polygon_t newpolygon;
+    newpolygon += Vector(1.,0.,0.);
+    newpolygon += Vector(0.,1.,0.);
+    newpolygon += Vector(0.,0.,1.);
+    SphericalPointDistribution::Polygon_t expectedpolygon;
+    expectedpolygon += Vector(1.,0.,0.);
+    expectedpolygon += Vector(0.,M_SQRT1_2,M_SQRT1_2);
+    SphericalPointDistribution::IndexTupleList_t matching;
+    SphericalPointDistribution::IndexList_t joined;
+    joined += 1,2;
+    matching += SphericalPointDistribution::IndexList_t(1,0);
+    matching += joined;
+    SphericalPointDistribution::IndexList_t IndexList =
+        SphericalPointDistribution::joinPoints(
+            newpolygon,
+            SphericalPointDistribution::VectorArray_t(newpolygon.begin(), newpolygon.end()),
+            matching);
+    SphericalPointDistribution::IndexList_t expected;
+    expected += 0,1;
+    CPPUNIT_ASSERT_EQUAL( expected, IndexList );
+    CPPUNIT_ASSERT_EQUAL( expectedpolygon, newpolygon );
+  }
+
+  // test with simple configuration of six points, two are joined, jumbled indices
+  {
+    SphericalPointDistribution::Polygon_t newpolygon;
+    newpolygon += Vector(1.,0.,1.);
+    newpolygon += Vector(1.,0.,0.);
+    newpolygon += Vector(1.,1.,0.);
+    newpolygon += Vector(0.,1.,0.);
+    newpolygon += Vector(0.,0.,1.);
+    newpolygon += Vector(1.,0.,1.);
+    SphericalPointDistribution::Polygon_t expectedpolygon;
+    expectedpolygon += Vector(1.,0.,1.);
+    expectedpolygon += Vector(1.,0.,0.);
+    expectedpolygon += Vector(1.,1.,0.);
+    expectedpolygon += Vector(1.,0.,1.);
+    expectedpolygon += Vector(0.,M_SQRT1_2,M_SQRT1_2); // new centers go last
+    SphericalPointDistribution::IndexTupleList_t matching;
+    SphericalPointDistribution::IndexList_t joined;
+    joined += 3,4;
+    matching += SphericalPointDistribution::IndexList_t(1,1);
+    matching += joined;
+    SphericalPointDistribution::IndexList_t IndexList =
+        SphericalPointDistribution::joinPoints(
+            newpolygon,
+            SphericalPointDistribution::VectorArray_t(newpolygon.begin(), newpolygon.end()),
+            matching);
+    SphericalPointDistribution::IndexList_t expected;
+    expected += 1,4;
+    CPPUNIT_ASSERT_EQUAL( expected, IndexList );
+    CPPUNIT_ASSERT_EQUAL( expectedpolygon, newpolygon );
+  }
+}
 
 /** UnitTest for matchSphericalPointDistributions() with three points

src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp

@@ -23 +23 @@
     CPPUNIT_TEST_SUITE( SphericalPointDistributionTest) ;
     CPPUNIT_TEST ( areEqualToWithinBoundsTest );
+    CPPUNIT_TEST ( joinPointsTest );
     CPPUNIT_TEST ( matchSphericalPointDistributionsTest_2 );
     CPPUNIT_TEST ( matchSphericalPointDistributionsTest_3 );
@@ -36 +37 @@
       void tearDown();
       void areEqualToWithinBoundsTest();
+      void joinPointsTest();
       void matchSphericalPointDistributionsTest_2();
       void matchSphericalPointDistributionsTest_3();