Changeset 0983e6

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:

3678eb

Parents:

653cea

git-author:

Frederik Heber <heber@…> (07/09/14 21:14:53)

git-committer:

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

Message:

Refactoring SphericalPointDistribution to combine point list and subset of indices.

• added struct PolygonWithIndices.

Location:

src/Fragmentation/Exporters

Files:

• SphericalPointDistribution.cpp (modified) (20 diffs)
• SphericalPointDistribution.hpp (modified) (3 diffs)

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

@@ -99 +99 @@
 }
 
+/** 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::PolygonWithIndices &_polygon)
+{
+  return calculateCenter(_polygon.polygon, _polygon.indices);
+}
+
 inline
 DistanceArray_t calculatePairwiseDistances(
@@ -177 +188 @@
     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);
+    const SphericalPointDistribution::PolygonWithIndices &_referencepositions,
+    const SphericalPointDistribution::PolygonWithIndices &_currentpositions)
+{
+  _oldCenter = calculateCenter(_referencepositions.polygon, _referencepositions.indices);
   // 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);
+  _newCenter = calculateCenter( _currentpositions.polygon, _currentpositions.indices);
 }
 /** Returns squared L2 error of the given \a _Matching.
@@ -251 +254 @@
 
 SphericalPointDistribution::Polygon_t SphericalPointDistribution::removeMatchingPoints(
-    const VectorArray_t &_points,
-    const IndexList_t &_matchingindices
+    const PolygonWithIndices &_points
     )
 {
-  SphericalPointDistribution::Polygon_t remainingreturnpolygon;
-  IndexArray_t indices(_matchingindices.begin(), _matchingindices.end());
+  SphericalPointDistribution::Polygon_t remainingpoints;
+  IndexArray_t indices(_points.indices.begin(), _points.indices.end());
   std::sort(indices.begin(), indices.end());
   LOG(4, "DEBUG: sorted matching is " << indices);
-  IndexArray_t remainingindices(_points.size(), -1);
+  IndexArray_t remainingindices(_points.polygon.size(), -1);
   std::generate(remainingindices.begin(), remainingindices.end(), UniqueNumber);
   IndexArray_t::iterator remainiter = std::set_difference(
@@ -269 +271 @@
   for (IndexArray_t::const_iterator iter = remainingindices.begin();
       iter != remainingindices.end(); ++iter) {
-    remainingreturnpolygon.push_back(_points[*iter]);
-  }
-
-  return remainingreturnpolygon;
+    remainingpoints.push_back(_points.polygon[*iter]);
+  }
+
+  return remainingpoints;
 }
 
@@ -501 +503 @@
 
 SphericalPointDistribution::Rotation_t SphericalPointDistribution::findPlaneAligningRotation(
-    const VectorArray_t &_referencepositions,
-    const VectorArray_t &_currentpositions,
-    const IndexList_t &_bestmatching
+    const PolygonWithIndices &_referencepositions,
+    const PolygonWithIndices &_currentpositions
     )
 {
@@ -520 +521 @@
   calculateOldAndNewCenters(
       oldCenter, newCenter,
-      _referencepositions, _currentpositions, _bestmatching);
+      _referencepositions, _currentpositions);
 
   if ((!oldCenter.IsZero()) && (!newCenter.IsZero())) {
@@ -538 +539 @@
     if ((oldCenter.IsZero()) && (newCenter.IsZero())) {
       // either oldCenter or newCenter (or both) is directly at origin
-      if (_bestmatching.size() == 2) {
+      if (_currentpositions.indices.size() == 2) {
         // line case
-        Vector oldPosition = _currentpositions[*_bestmatching.begin()];
-        Vector newPosition = _referencepositions[0];
+        Vector oldPosition = _currentpositions.polygon[*_currentpositions.indices.begin()];
+        Vector newPosition = _referencepositions.polygon[0];
         // check whether we need to rotate at all
         if (!oldPosition.IsEqualTo(newPosition)) {
@@ -555 +556 @@
         // both triangles/planes have same center, hence get axis by
         // VectorProduct of Normals
-        Plane newplane(_referencepositions[0], _referencepositions[1], _referencepositions[2]);
+        Plane newplane(_referencepositions.polygon[0], _referencepositions.polygon[1], _referencepositions.polygon[2]);
         VectorArray_t vectors;
-        for (IndexList_t::const_iterator iter = _bestmatching.begin();
-            iter != _bestmatching.end(); ++iter)
-          vectors.push_back(_currentpositions[*iter]);
+        for (IndexList_t::const_iterator iter = _currentpositions.indices.begin();
+            iter != _currentpositions.indices.end(); ++iter)
+          vectors.push_back(_currentpositions.polygon[*iter]);
         Plane oldplane(vectors[0], vectors[1], vectors[2]);
         Vector oldPosition = oldplane.getNormal();
@@ -607 +608 @@
 
 SphericalPointDistribution::Rotation_t SphericalPointDistribution::findPointAligningRotation(
-    const VectorArray_t &remainingold,
-    const VectorArray_t &remainingnew,
-    const IndexList_t &_bestmatching)
+    const PolygonWithIndices &remainingold,
+    const PolygonWithIndices &remainingnew)
 {
   // initialize rotation to zero
@@ -622 +622 @@
   calculateOldAndNewCenters(
       oldCenter, newCenter,
-      remainingold, remainingnew, _bestmatching);
-
-  Vector oldPosition = remainingnew[*_bestmatching.begin()];
-  Vector newPosition = remainingold[0];
-  LOG(6, "DEBUG: oldPosition is " << oldPosition << " and newPosition is " << newPosition);
+      remainingold, remainingnew);
+
+  Vector oldPosition = remainingnew.polygon[*remainingnew.indices.begin()];
+  Vector newPosition = remainingold.polygon[0];
+  LOG(6, "DEBUG: oldPosition is " << oldPosition << " (length: "
+      << oldPosition.Norm() << ") and newPosition is " << newPosition << " length(: "
+      << newPosition.Norm() << ")");
   if (!oldPosition.IsEqualTo(newPosition)) {
     if ((!oldCenter.IsZero()) && (!newCenter.IsZero())) {
-      oldCenter.Normalize();  // note weighted sum of normalized weight is not normalized
+      // we rotate at oldCenter and around the radial direction, which is again given
+      // by oldCenter.
       Rotation.first = oldCenter;
-      LOG(6, "DEBUG: Picking normalized oldCenter as Rotation.first " << oldCenter);
-      oldPosition.ProjectOntoPlane(Rotation.first);
-      newPosition.ProjectOntoPlane(Rotation.first);
+      Rotation.first.Normalize();  // note weighted sum of normalized weight is not normalized
+      LOG(6, "DEBUG: Using oldCenter " << oldCenter << " as rotation center and "
+          << Rotation.first << " as axis.");
+      oldPosition -= oldCenter;
+      newPosition -= oldCenter;
+      oldPosition = (oldPosition - oldPosition.Projection(Rotation.first));
+      newPosition = (newPosition - newPosition.Projection(Rotation.first));
       LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
     } else {
-      if (_bestmatching.size() == 2) {
+      if (remainingnew.indices.size() == 2) {
         // line situation
         try {
@@ -657 +664 @@
       } else {
         // triangle situation
-        Plane oldplane(remainingold[0], remainingold[1], remainingold[2]);
+        Plane oldplane(remainingold.polygon[0], remainingold.polygon[1], remainingold.polygon[2]);
         Rotation.first = oldplane.getNormal();
         LOG(6, "DEBUG: oldPlane is " << oldplane << " and normal is " << Rotation.first);
@@ -682 +689 @@
 {
   SphericalPointDistribution::Polygon_t remainingpoints;
-  VectorArray_t remainingold;
-  for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
-      iter != _polygon.end(); ++iter)
-    remainingold.push_back(iter->first);
+
   LOG(2, "INFO: Matching old polygon " << _polygon
       << " with new polygon " << _newpolygon);
@@ -699 +703 @@
     LOG(2, "INFO: Best matching is " << bestmatching);
 
-    // _newpolygon has changed, so now convert to array
-    VectorArray_t remainingnew(_newpolygon.begin(), _newpolygon.end());
+    const size_t NumberIds = std::min(bestmatching.size(), (size_t)3);
+    // create old set
+    PolygonWithIndices oldSet;
+    oldSet.indices.resize(NumberIds, -1);
+    std::generate(oldSet.indices.begin(), oldSet.indices.end(), UniqueNumber);
+    for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
+        iter != _polygon.end(); ++iter)
+      oldSet.polygon.push_back(iter->first);
+
+    // _newpolygon has changed, so now convert to array with matched indices
+    PolygonWithIndices newSet;
+    SphericalPointDistribution::IndexList_t::const_iterator beginiter = bestmatching.begin();
+    SphericalPointDistribution::IndexList_t::const_iterator enditer = bestmatching.begin();
+    std::advance(enditer, NumberIds);
+    newSet.indices.resize(NumberIds, -1);
+    std::copy(beginiter, enditer, newSet.indices.begin());
+    std::copy(_newpolygon.begin(),_newpolygon.end(), std::back_inserter(newSet.polygon));
 
     // determine rotation angles to align the two point distributions with
@@ -706 +725 @@
     // we use the center between the three first matching points
     /// the first rotation brings these two centers to coincide
-    VectorArray_t rotated_newpolygon = remainingnew;
+    PolygonWithIndices rotatednewSet = newSet;
     {
-      Rotation_t Rotation = findPlaneAligningRotation(
-          remainingold,
-          remainingnew,
-          bestmatching);
+      Rotation_t Rotation = findPlaneAligningRotation(oldSet, newSet);
       LOG(5, "DEBUG: Rotating coordinate system by " << Rotation.second
           << " around axis " << Rotation.first);
@@ -717 +733 @@
 
       // apply rotation angle to bring newCenter to oldCenter
-      for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
-          iter != rotated_newpolygon.end(); ++iter) {
+      for (VectorArray_t::iterator iter = rotatednewSet.polygon.begin();
+          iter != rotatednewSet.polygon.end(); ++iter) {
         Vector &current = *iter;
         LOG(6, "DEBUG: Original point is " << current);
@@ -732 +748 @@
         calculateOldAndNewCenters(
             oldCenter, rotatednewCenter,
-            remainingold, rotated_newpolygon, bestmatching);
+            oldSet, rotatednewSet);
         // NOTE: Center must not necessarily lie on the sphere with norm 1, hence, we
         // have to normalize it just as before, as oldCenter and newCenter lengths may differ.
@@ -740 +756 @@
           LOG(4, "CHECK: rotatednewCenter is " << rotatednewCenter
               << ", oldCenter is " << oldCenter);
-          ASSERT( (rotatednewCenter - oldCenter).NormSquared() < std::numeric_limits<double>::epsilon()*1e4,
+          const double difference = (rotatednewCenter - oldCenter).NormSquared();
+          ASSERT( difference < std::numeric_limits<double>::epsilon()*1e4,
               "matchSphericalPointDistributions() - rotation does not work as expected by "
-              +toString((rotatednewCenter - oldCenter).NormSquared())+".");
+              +toString(difference)+".");
         }
       }
@@ -750 +767 @@
     /// points themselves coincide
     if (bestmatching.size() > 1) {
-      Rotation_t Rotation = findPointAligningRotation(
-          remainingold,
-          rotated_newpolygon,
-          bestmatching);
+      Rotation_t Rotation = findPointAligningRotation(oldSet, rotatednewSet);
 
       // construct RotationAxis and two points on its plane, defining the angle
@@ -767 +781 @@
       {
         const IndexList_t::const_iterator iter = bestmatching.begin();
+
+        // check whether both old and newPosition are at same distance to oldCenter
+        Vector oldCenter = calculateCenter(oldSet);
+        const double distance = fabs(
+              (oldSet.polygon[0] - oldCenter).NormSquared()
+              - (rotatednewSet.polygon[*iter] - oldCenter).NormSquared()
+            );
+        LOG(4, "CHECK: Squared distance between oldPosition and newPosition "
+            << " with respect to oldCenter " << oldCenter << " is " << distance);
+//        ASSERT( distance < warn_amplitude,
+//            "matchSphericalPointDistributions() - old and newPosition's squared distance to oldCenter differs by "
+//            +toString(distance));
+
         Vector rotatednew = RotationAxis.rotateVector(
-            rotated_newpolygon[*iter],
+            rotatednewSet.polygon[*iter],
             Rotation.second);
         LOG(4, "CHECK: rotated first new bestmatching is " << rotatednew
-            << " while old was " << remainingold[0]);
-        ASSERT( (rotatednew - remainingold[0]).Norm() < warn_amplitude,
-            "matchSphericalPointDistributions() - orientation rotation ends up off by more than "
-            +toString(warn_amplitude)+".");
+            << " while old was " << oldSet.polygon[0]);
+        const double difference = (rotatednew - oldSet.polygon[0]).NormSquared();
+        ASSERT( difference < distance+1e-8,
+            "matchSphericalPointDistributions() - orientation rotation ends up off by "
+            +toString(difference)+", more than "
+            +toString(distance+1e-8)+".");
       }
 #endif
 
-      for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
-          iter != rotated_newpolygon.end(); ++iter) {
+      for (VectorArray_t::iterator iter = rotatednewSet.polygon.begin();
+          iter != rotatednewSet.polygon.end(); ++iter) {
         Vector &current = *iter;
         LOG(6, "DEBUG: Original point is " << current);
@@ -789 +818 @@
     // remove all points in matching and return remaining ones
     SphericalPointDistribution::Polygon_t remainingpoints =
-        removeMatchingPoints(rotated_newpolygon, bestmatching);
+        removeMatchingPoints(rotatednewSet);
     LOG(2, "INFO: Remaining points are " << remainingpoints);
     return remainingpoints;

src/Fragmentation/Exporters/SphericalPointDistribution.hpp

@@ -105 +105 @@
   static const double warn_amplitude;
 
+  struct PolygonWithIndices
+  {
+    //!> array with points
+    VectorArray_t polygon;
+    //!> list with indices for the above points, defining subset
+    IndexList_t indices;
+  };
+
 private:
   //!> grant unit tests access to private parts
@@ -115 +123 @@
 
   static Polygon_t removeMatchingPoints(
-      const VectorArray_t &_points,
-      const IndexList_t &_matchingindices
-      );
+      const PolygonWithIndices &_points);
 
   struct MatchingControlStructure {
@@ -149 +155 @@
 
   static Rotation_t findPlaneAligningRotation(
-      const VectorArray_t &_referencepositions,
-      const VectorArray_t &_currentpositions,
-      const IndexList_t &_bestmatching
+      const PolygonWithIndices &_referencepositions,
+      const PolygonWithIndices &_currentpositions
       );
 
   static Rotation_t findPointAligningRotation(
-      const VectorArray_t &remainingold,
-      const VectorArray_t &remainingnew,
-      const IndexList_t &_bestmatching);
+      const PolygonWithIndices &remainingold,
+      const PolygonWithIndices &remainingnew);
 
 };