Changeset d468b5 for src/Fragmentation/Exporters/SphericalPointDistribution.cpp

Timestamp:

May 25, 2016, 7:13:57 AM (9 years ago)

Author:

Frederik Heber <heber@…>

Children:

81557e

Parents:

3b6956

git-author:

Frederik Heber <heber@…> (05/10/16 18:33:50)

git-committer:

Frederik Heber <heber@…> (05/25/16 07:13:57)

Message:

Revert "Trimmed down SphericalPointDistribution to what is needed at the moment."

This reverts commit 0241c59e7f10c73f87336715af0aadac2b2d4db1.

• TESTFIX: SaturateAtoms regression tests are marked as XFAIL.
• commented-out addHydrogen section in SaturateAction.

File:

• src/Fragmentation/Exporters/SphericalPointDistribution.cpp (modified) (2 diffs)

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

@@ -38 +38 @@
 
 #include "CodePatterns/Assert.hpp"
-#include "CodePatterns/Log.hpp"
+#include "CodePatterns/IteratorAdaptors.hpp"
+#include "CodePatterns/toString.hpp"
 
 #include <algorithm>
+#include <cmath>
+#include <limits>
+#include <list>
 #include <vector>
+#include <map>
 
 #include "LinearAlgebra/Line.hpp"
@@ -47 +52 @@
 #include "LinearAlgebra/Vector.hpp"
 
-void SphericalPointDistribution::initSelf(const int _NumberOfPoints)
-{
-  switch (_NumberOfPoints)
-  {
-    case 0:
-      points = get<0>();
-      break;
-    case 1:
-      points = get<1>();
-      break;
-    case 2:
-      points = get<2>();
-      break;
-    case 3:
-      points = get<3>();
-      break;
-    case 4:
-      points = get<4>();
-      break;
-    case 5:
-      points = get<5>();
-      break;
-    case 6:
-      points = get<6>();
-      break;
-    case 7:
-      points = get<7>();
-      break;
-    case 8:
-      points = get<8>();
-      break;
-    case 9:
-      points = get<9>();
-      break;
-    case 10:
-      points = get<10>();
-      break;
-    case 11:
-      points = get<11>();
-      break;
-    case 12:
-      points = get<12>();
-      break;
-    case 14:
-      points = get<14>();
-      break;
-    default:
-      ASSERT(0, "SphericalPointDistribution::initSelf() - cannot deal with the case "
-          +toString(_NumberOfPoints)+".");
-  }
-  LOG(3, "DEBUG: Ideal polygon is " << points);
-}
-
+typedef std::list<unsigned int> IndexList_t;
+typedef std::vector<unsigned int> IndexArray_t;
+typedef std::vector<Vector> VectorArray_t;
+typedef std::vector<double> DistanceArray_t;
+
+DistanceArray_t calculatePairwiseDistances(
+    const std::vector<Vector> &_points,
+    const IndexList_t &_indices
+    )
+{
+  DistanceArray_t result;
+  for (IndexList_t::const_iterator firstiter = _indices.begin();
+      firstiter != _indices.end(); ++firstiter) {
+    for (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;
+}
+
+// class generator: taken from www.cplusplus.com example std::generate
+struct c_unique {
+  int current;
+  c_unique() {current=0;}
+  int operator()() {return ++current;}
+} UniqueNumber;
+
+/** Returns squared L2 error of the given \a _Matching.
+ *
+ * We compare the pair-wise distances of each associated matching
+ * and check whether these distances each match between \a _old and
+ * \a _new.
+ *
+ * \param _old first set of points (fewer or equal to \a _new)
+ * \param _new second set of points
+ * \param _Matching matching between the two sets
+ * \return pair with L1 and squared L2 error
+ */
+std::pair<double, double> calculateErrorOfMatching(
+    const std::vector<Vector> &_old,
+    const std::vector<Vector> &_new,
+    const IndexList_t &_Matching)
+{
+  std::pair<double, double> errors( std::make_pair( 0., 0. ) );
+
+  if (_Matching.size() > 1) {
+    // convert matching into two vectors to calculate distance among another
+
+    // calculate all pair-wise distances
+    IndexList_t keys(_Matching.size());
+    std::generate (keys.begin(), keys.end(), UniqueNumber);
+    const DistanceArray_t firstdistances = calculatePairwiseDistances(_old, keys);
+    const DistanceArray_t seconddistances = calculatePairwiseDistances(_new, _Matching);
+
+    ASSERT( firstdistances.size() == seconddistances.size(),
+        "calculateL2ErrorOfMatching() - mismatch in pair-wise distance array sizes.");
+    DistanceArray_t::const_iterator firstiter = firstdistances.begin();
+    DistanceArray_t::const_iterator seconditer = seconddistances.begin();
+    for (;(firstiter != firstdistances.end()) && (seconditer != seconddistances.end());
+        ++firstiter, ++seconditer) {
+      const double gap = *firstiter - *seconditer;
+      // L1 error
+      if (errors.first < gap)
+        errors.first = gap;
+      // L2 error
+      errors.second += gap*gap;
+    }
+  }
+
+  return errors;
+}
+
+SphericalPointDistribution::Polygon_t removeMatchingPoints(
+    const SphericalPointDistribution::Polygon_t &_points,
+    const IndexList_t &_matchingindices
+    )
+{
+  SphericalPointDistribution::Polygon_t remainingpoints;
+  IndexArray_t indices(_matchingindices.begin(), _matchingindices.end());
+  std::sort(indices.begin(), indices.end());
+  IndexArray_t::const_iterator valueiter = indices.begin();
+  SphericalPointDistribution::Polygon_t::const_iterator pointiter =
+      _points.begin();
+  for (unsigned int i=0; i< _points.size(); ++i, ++pointiter) {
+    // skip all those in values
+    if (*valueiter == i)
+      ++valueiter;
+    else
+      remainingpoints.push_back(*pointiter);
+  }
+
+  return remainingpoints;
+}
+
+/** Rotates a given polygon around x, y, and z axis by the given angles.
+ *
+ * Essentially, we concentrate on the three points of the polygon to rotate
+ * to the correct position. First, we rotate its center via \a angles,
+ * then we rotate the "triangle" around itself/\a _RotationAxis by
+ * \a _RotationAngle.
+ *
+ * \param _polygon polygon whose points to rotate
+ * \param _angles vector with rotation angles for x,y,z axis
+ * \param _RotationAxis
+ * \param _RotationAngle
+ */
+SphericalPointDistribution::Polygon_t rotatePolygon(
+    const SphericalPointDistribution::Polygon_t &_polygon,
+    const std::vector<double> &_angles,
+    const Line &_RotationAxis,
+    const double _RotationAngle)
+{
+  SphericalPointDistribution::Polygon_t rotated_polygon = _polygon;
+  RealSpaceMatrix rotation;
+  ASSERT( _angles.size() == 3,
+      "rotatePolygon() - not exactly "+toString(3)+" angles given.");
+  rotation.setRotation(_angles[0] * M_PI/180., _angles[1] * M_PI/180., _angles[2] * M_PI/180.);
+
+  // apply rotation angles
+  for (SphericalPointDistribution::Polygon_t::iterator iter = rotated_polygon.begin();
+      iter != rotated_polygon.end(); ++iter) {
+    *iter = rotation * (*iter);
+    _RotationAxis.rotateVector(*iter, _RotationAngle);
+  }
+
+  return rotated_polygon;
+}
+
+struct MatchingControlStructure {
+  bool foundflag;
+  double bestL2;
+  IndexList_t bestmatching;
+  VectorArray_t oldpoints;
+  VectorArray_t newpoints;
+};
+
+/** Recursive function to go through all possible matchings.
+ *
+ * \param _MCS structure holding global information to the recursion
+ * \param _matching current matching being build up
+ * \param _indices contains still available indices
+ * \param _matchingsize
+ */
+void recurseMatchings(
+    MatchingControlStructure &_MCS,
+    IndexList_t &_matching,
+    IndexList_t _indices,
+    unsigned int _matchingsize)
+{
+  //!> threshold for L1 error below which matching is immediately acceptable
+  const double L1THRESHOLD = 1e-2;
+  if (!_MCS.foundflag) {
+    if (_matching.size() < _matchingsize) {
+      // go through all indices
+      for (IndexList_t::iterator iter = _indices.begin();
+          iter != _indices.end();) {
+        // add index to matching
+        _matching.push_back(*iter);
+        // 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);
+        // re-add chosen index and reset index to new position
+        _indices.insert(backupiter, _matching.back());
+        iter = backupiter;
+        // remove index from _matching to make space for the next one
+        _matching.pop_back();
+      }
+      // gone through all indices then exit recursion
+      _MCS.foundflag = true;
+    } else {
+      // calculate errors
+      std::pair<double, double> errors = calculateErrorOfMatching(
+          _MCS.oldpoints, _MCS.newpoints, _matching);
+      if (errors.first < L1THRESHOLD) {
+        _MCS.bestmatching = _matching;
+        _MCS.foundflag = true;
+      }
+      if (_MCS.bestL2 > errors.second) {
+        _MCS.bestmatching = _matching;
+        _MCS.bestL2 = errors.second;
+      }
+    }
+  }
+}
+
+SphericalPointDistribution::Polygon_t
+SphericalPointDistribution::matchSphericalPointDistributions(
+    const SphericalPointDistribution::Polygon_t &_polygon,
+    const SphericalPointDistribution::Polygon_t &_newpolygon
+    )
+{
+  SphericalPointDistribution::Polygon_t remainingpoints;
+  VectorArray_t remainingold(_polygon.begin(), _polygon.end());
+  VectorArray_t remainingnew(_newpolygon.begin(), _newpolygon.end());
+
+  if (_polygon.size() > 0) {
+    MatchingControlStructure MCS;
+    MCS.foundflag = false;
+    MCS.bestL2 = std::numeric_limits<double>::max();
+    MCS.oldpoints.insert(MCS.oldpoints.begin(), _polygon.begin(),_polygon.end() );
+    MCS.newpoints.insert(MCS.newpoints.begin(), _newpolygon.begin(),_newpolygon.end() );
+
+    // search for bestmatching combinatorially
+    {
+      // translate polygon into vector to enable index addressing
+      IndexList_t indices(_newpolygon.size());
+      std::generate(indices.begin(), indices.end(), UniqueNumber);
+      IndexList_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);
+    }
+
+    // determine rotation angles to align the two point distributions with
+    // respect to bestmatching
+    std::vector<double> angles(3);
+    Vector newCenter;
+    {
+      // calculate center of triangle/line/point consisting of first points of matching
+      Vector oldCenter;
+      IndexList_t::const_iterator iter = MCS.bestmatching.begin();
+      unsigned int i = 0;
+      for (; (i<3) && (i<MCS.bestmatching.size()); ++i, ++iter) {
+        oldCenter += remainingold[i];
+        newCenter += remainingnew[*iter];
+      }
+      oldCenter *= 1./(double)i;
+      newCenter *= 1./(double)i;
+
+      Vector direction(0.,0.,0.);
+      for(size_t i=0;i<NDIM;++i) {
+        // create new rotation axis
+        direction[i] = 1.;
+        const Line axis (zeroVec, direction);
+        // calculate rotation angle for this axis
+        const double alpha = direction.Angle(oldCenter) - direction.Angle(newCenter);
+        // perform rotation
+        axis.rotateVector(newCenter, alpha);
+        // store angle
+        angles[i] = alpha;
+        // reset direction component for next iteration
+        direction[i] = 0.;
+      }
+    }
+    const Line RotationAxis(zeroVec, newCenter);
+    const double RotationAngle =
+        newCenter.Angle(remainingold[0])
+        - newCenter.Angle(remainingnew[*MCS.bestmatching.begin()]);
+
+    // rotate _newpolygon
+    SphericalPointDistribution::Polygon_t rotated_newpolygon =
+        rotatePolygon(_newpolygon, angles, RotationAxis, RotationAngle);
+
+    // remove all points in matching and return remaining ones
+    return removeMatchingPoints(rotated_newpolygon, MCS.bestmatching);
+  } else
+    return _newpolygon;
+}
+
+
+