Changeset 4f04ab8 for src/Graph/BondGraph.cpp

Timestamp:

Apr 18, 2013, 4:23:26 PM (12 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, 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_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

53c579

Parents:

326000

git-author:

Frederik Heber <heber@…> (03/21/13 09:33:31)

git-committer:

Frederik Heber <heber@…> (04/18/13 16:23:26)

Message:

Using boost::graph to obtain maximum matching for bond degree correction.

• added bimap to boost.m4, added boost_bimap and boost_graph to configure.ac.
• rewrote createMatching() in BondGraph with boost::graph, code mostly taken from maximum cardinality matching example.

File:

• src/Graph/BondGraph.cpp (modified) (2 diffs)

src/Graph/BondGraph.cpp

@@ -33 +33 @@
 #endif
 
+// boost::graph uses placement new
+#include <boost/graph/adjacency_list.hpp>
+
 #include "CodePatterns/MemDebug.hpp"
 
+#include <algorithm>
+#include <boost/bimap.hpp>
+#include <boost/bind.hpp>
+#include <boost/foreach.hpp>
+#include <boost/function.hpp>
+#include <boost/graph/max_cardinality_matching.hpp>
 #include <deque>
 #include <iostream>
@@ -320 +329 @@
   return treeedges;
 }
+
+struct hasDelta {
+  bool operator()(atom *_atom) {
+    const double delta =
+        _atom->getType()->getNoValenceOrbitals() - _atom->CountBonds();
+    return (fabs(delta) > 0);
+  }
+};
+
+typedef std::set<atom *> deltaatoms_t;
+typedef std::set<bond::ptr> deltaedges_t;
+
+static int gatherAllDeltaAtoms(
+    const deltaatoms_t &allatoms,
+    deltaatoms_t &deltaatoms)
+{
+  static hasDelta delta;
+  deltaatoms.clear();
+  for (deltaatoms_t::const_iterator iter = allatoms.begin();
+      iter != allatoms.end(); ++iter) {
+    if (delta(*iter))
+      deltaatoms.insert(*iter);
+  }
+  return deltaatoms.size();
+}
+
+typedef boost::bimap<int,atom*> AtomIndexLookup_t;
+
+static AtomIndexLookup_t createAtomIndexLookup(
+    const deltaedges_t &edges)
+{
+  AtomIndexLookup_t AtomIndexLookup;
+  size_t index = 0;
+  for (deltaedges_t::const_iterator edgeiter = edges.begin();
+      edgeiter != edges.end(); ++edgeiter) {
+    const bond::ptr & _bond = *edgeiter;
+
+    // insert left into lookup
+    std::pair< AtomIndexLookup_t::iterator, bool > inserter =
+        AtomIndexLookup.insert( AtomIndexLookup_t::value_type(index, _bond->leftatom) );
+    if (inserter.second)
+      ++index;
+
+    // insert right into lookup
+    inserter = AtomIndexLookup.insert( AtomIndexLookup_t::value_type(index, _bond->rightatom) );
+    if (inserter.second)
+      ++index;
+  }
+  return AtomIndexLookup;
+}
+
+typedef std::map< std::pair<atom *, atom*>, bond::ptr> BondLookup_t;
+static BondLookup_t createBondLookup(
+    const deltaedges_t &edges)
+{
+  BondLookup_t BondLookup;
+  for (deltaedges_t::const_iterator edgeiter = edges.begin();
+      edgeiter != edges.end(); ++edgeiter) {
+    const bond::ptr & _bond = *edgeiter;
+
+    // insert bond into pair lookup
+    BondLookup.insert(
+        std::make_pair(
+            std::make_pair(_bond->leftatom, _bond->rightatom),
+            _bond)
+    );
+  }
+  return BondLookup;
+}
+
+typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> graph_t;
+typedef std::vector<boost::graph_traits<graph_t>::vertex_descriptor> mate_t;
+
+static void fillEdgesIntoMatching(
+    graph_t &g,
+    mate_t &mate,
+    const AtomIndexLookup_t &AtomIndexLookup,
+    const BondLookup_t &BondLookup,
+    deltaedges_t &matching
+    )
+{
+  matching.clear();
+  boost::graph_traits<graph_t>::vertex_iterator vi, vi_end;
+  for(tie(vi,vi_end) = boost::vertices(g); vi != vi_end; ++vi)
+    if (mate[*vi] != boost::graph_traits<graph_t>::null_vertex() && *vi < mate[*vi]) {
+      atom * leftatom = AtomIndexLookup.left.at(*vi);
+      atom * rightatom = AtomIndexLookup.left.at(mate[*vi]);
+      std::pair<atom *,atom *> AtomPair(leftatom, rightatom);
+      std::pair<atom *,atom *> reverseAtomPair(rightatom, leftatom);
+      BondLookup_t::const_iterator iter = BondLookup.find(AtomPair);
+      if (iter != BondLookup.end()) {
+        matching.insert(iter->second);
+      } else {
+        iter = BondLookup.find(reverseAtomPair);
+        if (iter != BondLookup.end()) {
+          matching.insert(iter->second);
+        } else
+          ASSERT(0, "fillEdgesIntoMatching() - cannot find ("+toString(*vi)+","
+              +toString(mate[*vi])+") in BondLookup.");
+      }
+    }
+}
+
+static bool createMatching(
+    deltaedges_t &deltaedges,
+    deltaedges_t &matching
+    )
+{
+  // gather all vertices for graph in a lookup structure
+  // to translate the graphs indices to atoms and also to associated bonds
+  AtomIndexLookup_t AtomIndexLookup = createAtomIndexLookup(deltaedges);
+  BondLookup_t BondLookup = createBondLookup(deltaedges);
+  const int n_vertices = AtomIndexLookup.size();
+
+  // construct graph
+  graph_t g(n_vertices);
+
+  // add edges to graph
+  for (deltaedges_t::const_iterator edgeiter = deltaedges.begin();
+      edgeiter != deltaedges.end(); ++edgeiter) {
+    const bond::ptr & _bond = *edgeiter;
+    boost::add_edge(
+        AtomIndexLookup.right.at(_bond->leftatom),
+        AtomIndexLookup.right.at(_bond->rightatom),
+        g);
+  }
+
+  // mate structure contains unique edge partner to every vertex in matching
+  mate_t mate(n_vertices);
+
+  // get maximum matching over given edges
+  bool success = boost::checked_edmonds_maximum_cardinality_matching(g, &mate[0]);
+
+  if (success) {
+    LOG(1, "STATUS: Found a matching of size " << matching_size(g, &mate[0]) << ".");
+    fillEdgesIntoMatching(g, mate, AtomIndexLookup, BondLookup, matching);
+  } else {
+    ELOG(2, "Failed to find maximum matching.");
+  }
+
+  return success;
+}
+
+
+int BondGraph::calculateBondDegree(const deltaatoms_t &allatoms) const
+{
+  deltaatoms_t deltaatoms;
+  deltaedges_t deltaedges;
+  deltaedges_t matching;
+
+  LOG(1, "STATUS: Calculating bond degrees.");
+  if (DoLog(2)) {
+    std::stringstream output;
+    output << "INFO: All atoms are: ";
+    BOOST_FOREACH( atom *_atom, allatoms) {
+      output << *_atom << " ";
+    }
+    LOG(2, output.str());
+  }
+
+  size_t IterCount = 0;
+  // 1. gather all atoms with delta > 0
+  while ((gatherAllDeltaAtoms(allatoms, deltaatoms) != 0) && (IterCount < 100)) {
+    if (DoLog(3)) {
+      std::stringstream output;
+      output << "DEBUG: Current delta atoms are: ";
+      BOOST_FOREACH( atom *_atom, deltaatoms) {
+        output << *_atom << " ";
+      }
+      LOG(3, output.str());
+    }
+
+    // 2. gather all edges that connect atoms with both(!) having delta > 0
+    deltaedges.clear();
+    for (deltaatoms_t::const_iterator atomiter = deltaatoms.begin();
+        atomiter != deltaatoms.end(); ++atomiter) {
+      const atom * const _atom = *atomiter;
+      const BondList& ListOfBonds = (*atomiter)->getListOfBonds();
+      for (BondList::const_iterator bonditer = ListOfBonds.begin();
+          bonditer != ListOfBonds.end();++bonditer) {
+        const bond::ptr _bond = *bonditer;
+        if ((_bond->leftatom == _atom) && (deltaatoms.find(_bond->rightatom) != deltaatoms.end()))
+          deltaedges.insert(_bond);
+        else if ((_bond->rightatom == _atom) && (deltaatoms.find(_bond->leftatom) != deltaatoms.end()))
+          deltaedges.insert(_bond);
+      }
+    }
+    if (DoLog(3)) {
+      std::stringstream output;
+      output << "DEBUG: Current delta bonds are: ";
+      BOOST_FOREACH( bond::ptr _bond, deltaedges) {
+        output << *_bond << " ";
+      }
+      LOG(3, output.str());
+    }
+    if (deltaedges.empty())
+      break;
+
+    // 3. build matching over these edges
+    if (!createMatching(deltaedges, matching))
+      break;
+    if (DoLog(2)) {
+      std::stringstream output;
+      output << "DEBUG: Resulting matching is: ";
+      BOOST_FOREACH( bond::ptr _bond, matching) {
+        output << *_bond << " ";
+      }
+      LOG(2, output.str());
+    }
+
+    // 4. increase degree for each edge of the matching
+    LOG(2, "DEBUG: Increasing bond degrees by one.");
+    for (deltaedges_t::iterator edgeiter = matching.begin();
+        edgeiter != matching.end(); ++edgeiter) {
+      (*edgeiter)->setDegree( (*edgeiter)->getDegree()+1 );
+    }
+
+    // 6. stop after a given number of iterations or when all atoms are happy.
+    ++IterCount;
+  };
+
+  // check whether we still have deltaatoms
+  {
+    hasDelta delta;
+    for (deltaatoms_t::const_iterator iter = allatoms.begin();
+        iter != allatoms.end(); ++iter)
+      if (delta(*iter))
+        ELOG(2, "Could not find satisfy charge neutrality for atom " << **iter << ".");
+  }
+
+  return deltaedges.size();
+}