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

Timestamp:

Sep 14, 2016, 7:02:33 AM (9 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:

1af2ae

Parents:

36bd59

git-author:

Frederik Heber <heber@…> (07/12/14 19:07:44)

git-committer:

Frederik Heber <heber@…> (09/14/16 07:02:33)

Message:

recurseMatching() now takes connections into account.

• matchSphericalPointDistribution() replaced by getRemainingPoints() as that describes what it does. match...() sounds symmetrical which the function is no longer as connection is associated with former _newpolygon.
• SphericalPointDistribution now has internal points and adjacency, initialized by initSelf().
• findBestMatching() and getRemainingPoints() are no longer static functions.
• all unit tests are working again, including the .._multiple() that tests for joint points due to bond degree greater than 1.
• the huge case switch in SaturatedFragment::saturateAtom() now resides in initSelf().

File:

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

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

@@ -480 +480 @@
       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) {
@@ -490 +491 @@
           }
         }
-        // advance in matching to same position
+
+        // advance in matching to current bin to fill in
         const size_t OldIndex = std::distance(_remainingweights.begin(), _remainiter);
         while (_matching.size() <= OldIndex) { // add empty lists of new bin is opened
@@ -498 +500 @@
         IndexTupleList_t::iterator filliniter = _matching.begin();
         std::advance(filliniter, OldIndex);
+
+        // check whether connection between bins' indices and candidate is satisfied
+        {
+          adjacency_t::const_iterator finder = _MCS.adjacency.find(*iter);
+          ASSERT( finder != _MCS.adjacency.end(),
+              "recurseMatchings() - "+toString(*iter)+" is not in adjacency list.");
+          if ((!filliniter->empty())
+              && (finder->second.find(*filliniter->begin()) == finder->second.end())) {
+            LOG(5, "DEBUG; Skipping index " << *iter
+                << " as is not connected to current set." << *filliniter << ".");
+            ++iter; // note that for loop does not contain incrementor
+            continue;
+          }
+        }
+
         // add index to matching
         filliniter->push_back(*iter);
@@ -532 +549 @@
 }
 
+SphericalPointDistribution::MatchingControlStructure::MatchingControlStructure(
+    const adjacency_t &_adjacency,
+    const VectorArray_t &_oldpoints,
+    const VectorArray_t &_newpoints,
+    const WeightsArray_t &_weights
+    ) :
+  foundflag(false),
+  bestL2(std::numeric_limits<double>::max()),
+  adjacency(_adjacency),
+  oldpoints(_oldpoints),
+  newpoints(_newpoints),
+  weights(_weights)
+{}
+
 /** Finds combinatorially the best matching between points in \a _polygon
  * and \a _newpolygon.
@@ -551 +582 @@
  */
 SphericalPointDistribution::IndexList_t SphericalPointDistribution::findBestMatching(
-    const WeightedPolygon_t &_polygon,
-    Polygon_t &_newpolygon
+    const WeightedPolygon_t &_polygon
     )
 {
-  MatchingControlStructure MCS;
-  MCS.foundflag = false;
-  MCS.bestL2 = std::numeric_limits<double>::max();
   // transform lists into arrays
+  VectorArray_t oldpoints;
+  VectorArray_t newpoints;
+  WeightsArray_t weights;
   for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
       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() );
+    oldpoints.push_back(iter->first);
+    weights.push_back(iter->second);
+  }
+  newpoints.insert(newpoints.begin(), points.begin(), points.end() );
+  MatchingControlStructure MCS(adjacency, oldpoints, newpoints, weights);
 
   // search for bestmatching combinatorially
   {
     // translate polygon into vector to enable index addressing
-    IndexList_t indices(_newpolygon.size());
+    IndexList_t indices(points.size());
     std::generate(indices.begin(), indices.end(), UniqueNumber);
     IndexTupleList_t matching;
@@ -594 +625 @@
   // combine multiple points and create simple IndexList from IndexTupleList
   const SphericalPointDistribution::IndexList_t IndexList =
-      joinPoints(_newpolygon, MCS.newpoints, MCS.bestmatching);
+      joinPoints(points, MCS.newpoints, MCS.bestmatching);
 
   return IndexList;
@@ -781 +812 @@
 }
 
+void SphericalPointDistribution::initSelf(const int _NumberOfPoints)
+{
+  switch (_NumberOfPoints)
+  {
+    case 0:
+      points = get<0>();
+      adjacency = getConnections<0>();
+      break;
+    case 1:
+      points = get<1>();
+      adjacency = getConnections<1>();
+      break;
+    case 2:
+      points = get<2>();
+      adjacency = getConnections<2>();
+      break;
+    case 3:
+      points = get<3>();
+      adjacency = getConnections<3>();
+      break;
+    case 4:
+      points = get<4>();
+      adjacency = getConnections<4>();
+      break;
+    case 5:
+      points = get<5>();
+      adjacency = getConnections<5>();
+      break;
+    case 6:
+      points = get<6>();
+      adjacency = getConnections<6>();
+      break;
+    case 7:
+      points = get<7>();
+      adjacency = getConnections<7>();
+      break;
+    case 8:
+      points = get<8>();
+      adjacency = getConnections<8>();
+      break;
+    case 9:
+      points = get<9>();
+      adjacency = getConnections<9>();
+      break;
+    case 10:
+      points = get<10>();
+      adjacency = getConnections<10>();
+      break;
+    case 11:
+      points = get<11>();
+      adjacency = getConnections<11>();
+      break;
+    case 12:
+      points = get<12>();
+      adjacency = getConnections<12>();
+      break;
+    case 14:
+      points = get<14>();
+      adjacency = getConnections<14>();
+      break;
+    default:
+      ASSERT(0, "SphericalPointDistribution::initSelf() - cannot deal with the case "
+          +toString(_NumberOfPoints)+".");
+  }
+  LOG(3, "DEBUG: Ideal polygon is " << points);
+}
 
 SphericalPointDistribution::Polygon_t
-SphericalPointDistribution::matchSphericalPointDistributions(
-    const SphericalPointDistribution::WeightedPolygon_t &_polygon,
-    SphericalPointDistribution::Polygon_t &_newpolygon
-    )
+SphericalPointDistribution::getRemainingPoints(
+    const WeightedPolygon_t &_polygon,
+    const int _N)
 {
   SphericalPointDistribution::Polygon_t remainingpoints;
 
+  // initialze to given number of points
+  initSelf(_N);
   LOG(2, "INFO: Matching old polygon " << _polygon
-      << " with new polygon " << _newpolygon);
-
-  if (_polygon.size() == _newpolygon.size()) {
-    // same number of points desired as are present? Do nothing
-    LOG(2, "INFO: There are no vacant points to return.");
+      << " with new polygon " << points);
+
+  // check whether any points will remain vacant
+  int RemainingPoints = _N;
+  for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
+      iter != _polygon.end(); ++iter)
+    RemainingPoints -= iter->second;
+  if (RemainingPoints == 0)
     return remainingpoints;
-  }
-
-  if (_polygon.size() > 0) {
-    IndexList_t bestmatching = findBestMatching(_polygon, _newpolygon);
+
+  if (_N > 0) {
+    IndexList_t bestmatching = findBestMatching(_polygon);
     LOG(2, "INFO: Best matching is " << bestmatching);
 
@@ -819 +919 @@
     newSet.indices.resize(NumberIds, -1);
     std::copy(beginiter, enditer, newSet.indices.begin());
-    std::copy(_newpolygon.begin(),_newpolygon.end(), std::back_inserter(newSet.polygon));
+    std::copy(points.begin(),points.end(), std::back_inserter(newSet.polygon));
 
     // determine rotation angles to align the two point distributions with
@@ -927 +1027 @@
     return remainingpoints;
   } else
-    return _newpolygon;
-}
+    return points;
+}