Changeset 64cafb2 for src

Timestamp:

Sep 10, 2016, 4:01:20 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, Fixes, 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, SaturateAtoms_singleDegree, 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:

0c42f2, cdac1d

Parents:

d32f60

git-author:

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

git-committer:

Frederik Heber <heber@…> (09/10/16 16:01:20)

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.

Location:

src

Files:

• Actions/AtomAction/SaturateAction.cpp (modified) (1 diff)
• Fragmentation/Exporters/SphericalPointDistribution.cpp (modified) (2 diffs)
• Fragmentation/Exporters/SphericalPointDistribution.hpp (modified) (2 diffs)
• Fragmentation/Exporters/unittests/Makefile.am (modified) (3 diffs)
• Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp (added)
• Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp (added)

src/Actions/AtomAction/SaturateAction.cpp

@@ -77 +77 @@
         "ChangeElementAction::performCall() - mol is still NULL." );
 
-    // add the hydrogens
-    const Vector AtomsPosition = _atom->getPosition();
-    double typical_distance = _atom->getType()->getHBondDistance(0);
-    if (typical_distance == -1.)
-      typical_distance = 1.;
-    SphericalPointDistribution PointSphere(typical_distance);
-    PointSphere.initSelf(_atom->getType()->getNoValenceOrbitals());
-    for (SphericalPointDistribution::Polygon_t::const_iterator iter = PointSphere.points.begin();
-        iter != PointSphere.points.end(); ++iter) {
-      // for every Vector add a point relative to atom's position
-      atom * hydrogenAtom = World::getInstance().createAtom();
-      hydrogenAtom->setType(hydrogen);
-      hydrogenAtom->setPosition(AtomsPosition + *iter);
-      mol->AddAtom(hydrogenAtom);
-      // add bond
-      _atom->addBond(WorldTime::getTime(), hydrogenAtom);
-      // mark down for undo state
-      addedHydrogens.push_back(AtomicInfo(*(hydrogenAtom)));
-    }
+//    // add the hydrogens
+//    const Vector AtomsPosition = _atom->getPosition();
+//    double typical_distance = _atom->getType()->getHBondDistance(0);
+//    if (typical_distance == -1.)
+//      typical_distance = 1.;
+//    SphericalPointDistribution PointSphere(typical_distance);
+//    PointSphere.initSelf(_atom->getType()->getNoValenceOrbitals());
+//    for (SphericalPointDistribution::Polygon_t::const_iterator iter = PointSphere.points.begin();
+//        iter != PointSphere.points.end(); ++iter) {
+//      // for every Vector add a point relative to atom's position
+//      atom * hydrogenAtom = World::getInstance().createAtom();
+//      hydrogenAtom->setType(hydrogen);
+//      hydrogenAtom->setPosition(AtomsPosition + *iter);
+//      mol->AddAtom(hydrogenAtom);
+//      // add bond
+//      _atom->addBond(WorldTime::getTime(), hydrogenAtom);
+//      // mark down for undo state
+//      addedHydrogens.push_back(AtomicInfo(*(hydrogenAtom)));
+//    }
   }
 

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;
+}
+
+
+

src/Fragmentation/Exporters/SphericalPointDistribution.hpp

@@ -57 +57 @@
   }
 
-  /** Initializes the polygon with the given \a _NumberOfPoints.
+
+  /** Matches a given spherical distribution with another containing more
+   * points.
    *
-   * \param _NumberOfPoints number of points
+   * This is a helper to determine points where to best insert saturation
+   * hydrogens.
+   *
+   * \param _polygon current occupied positions
+   * \param _newpolygon ideal distribution to match best with current occupied
+   *        positions
+   * \return remaining vacant positions relative to \a _polygon
    */
-  void initSelf(const int _NumberOfPoints);
+  static Polygon_t matchSphericalPointDistributions(
+      const Polygon_t &_polygon,
+      const Polygon_t &_newpolygon
+      );
+
 
   //!> default radius of the spherical distribution
@@ -67 +79 @@
   //!> precalculated value for root of 3
   const double SQRT_3;
-  //!> points in this polygon
-  Polygon_t points;
 };
 

src/Fragmentation/Exporters/unittests/Makefile.am

@@ -5 +5 @@
         ../Fragmentation/Exporters/unittests/HydrogenPoolUnitTest.cpp \
         ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.cpp \
-        ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.cpp
+        ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.cpp \
+        ../Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp
 
 FRAGMENTATIONEXPORTERSTESTSHEADERS = \
         ../Fragmentation/Exporters/unittests/HydrogenPoolUnitTest.hpp \
         ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.hpp \
-        ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.hpp
+        ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.hpp \
+        ../Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp
 
 FRAGMENTATIONEXPORTERSTESTS = \
         HydrogenPoolUnitTest \
         SaturatedFragmentUnitTest \
-        SaturationDistanceMaximizerUnitTest
+        SaturationDistanceMaximizerUnitTest \
+        SphericalPointDistributionUnitTest
 
 TESTS += $(FRAGMENTATIONEXPORTERSTESTS)
@@ -45 +48 @@
 
 SaturationDistanceMaximizerUnitTest_SOURCES = \
+        ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.cpp \
         ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.hpp \
         ../Fragmentation/Exporters/unittests/stubs/SaturatedBondStub.cpp
@@ -51 +55 @@
         $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
         ${FRAGMENTATIONEXPORTERSLIBS}
+        
+SphericalPointDistributionUnitTest_SOURCES = \
+        ../Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp \
+        ../Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp
+SphericalPointDistributionUnitTest_LDADD = \
+        ../libMolecuilderUI.la  \
+        $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
+        ${FRAGMENTATIONEXPORTERSLIBS}