Changeset 6e5907

Timestamp:

Jul 5, 2017, 7:40:48 PM (8 years ago)

Author:

Frederik Heber <frederik.heber@…>

Branches:

ForceAnnealing_oldresults, IndependentFragmentGrids_IntegrationTest

Children:

89235ea

Parents:

6afd46

git-author:

Frederik Heber <frederik.heber@…> (05/18/17 17:45:47)

git-committer:

Frederik Heber <frederik.heber@…> (07/05/17 19:40:48)

Message:

Extracted extraction (subset of) nodes from BoostGraph into BreadthFirstSearchGatherer.

• also added helper namespace BoostGraphHelpers.
• we now treat the vertex indices and vertex names properly. Before that they had to coincide. Now, the name is the atomic id associated with the node and the index is the boost::graph internal index.

Location:

src

Files:

• Actions/MoleculeAction/StretchBondAction.cpp (modified) (8 diffs)
• Actions/MoleculeAction/StretchBondAction.def (modified) (1 diff)
• Graph/BoostGraphCreator.cpp (modified) (3 diffs)
• Graph/BoostGraphCreator.hpp (modified) (4 diffs)
• Graph/BoostGraphHelpers.hpp (added)
• Graph/BreadthFirstSearchGatherer.cpp (added)
• Graph/BreadthFirstSearchGatherer.hpp (added)
• Graph/Makefile.am (modified) (2 diffs)

src/Actions/MoleculeAction/StretchBondAction.cpp

@@ -33 +33 @@
 #endif
 
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/breadth_first_search.hpp>
-
 //#include "CodePatterns/MemDebug.hpp"
 
@@ -52 +49 @@
 #include "Descriptors/AtomIdDescriptor.hpp"
 #include "Graph/BoostGraphCreator.hpp"
+#include "Graph/BoostGraphHelpers.hpp"
+#include "Graph/BreadthFirstSearchGatherer.hpp"
 #include "molecule.hpp"
 #include "World.hpp"
@@ -61 +60 @@
 #include "Action_impl_pre.hpp"
 
-
-/**
- * I have no idea why this is so complicated with BGL ...
- *
- * This is taken from the book "The Boost Graph Library: User Guide and Reference Manual, Portable Documents",
- * chapter "Basic Graph Algorithms", example on calculating the bacon number.
- */
-template <typename DistanceMap>
-class distance_recorder : public boost::default_bfs_visitor
-{
-public:
-  distance_recorder(DistanceMap dist) : d(dist) {}
-
-  template <typename Edge, typename Graph>
-  void tree_edge(Edge e, const Graph &g) const {
-    typename boost::graph_traits<Graph>::vertex_descriptor u = source(e,g), v = target(e,g);
-    d[v] = d[u] + 1;
-  }
-
-private:
-  DistanceMap d;
-};
-
-template <typename DistanceMap>
-distance_recorder<DistanceMap> record_distance(DistanceMap d)
-{
-  return distance_recorder<DistanceMap>(d);
-}
 
 static bool addEdgePredicate(
@@ -110 +81 @@
     return Action::failure;
   }
+  molecule *mol = World::getInstance().
+      getMolecule(MoleculeById(atoms[0]->getMolecule()->getId()));
+  if (mol != atoms[1]->getMolecule()) {
+    STATUS("The two selected atoms must belong to the same molecule.");
+    return Action::failure;
+  }
+
+  // gather undo information
+  const double olddistance = atoms[0]->getPosition().distance(atoms[1]->getPosition());
+  const double newdistance = params.bonddistance.get();
+  LOG(1, "INFO: Old bond distance is " << olddistance << ", stretching to " << newdistance << ".");
+
+  // gather sorted ids
   std::vector<atomId_t> atomids(2);
   atomids[0] = atoms[0]->getId();
@@ -115 +99 @@
   std::sort(atomids.begin(), atomids.end());
   LOG(1, "DEBUG: Selected nodes are " << atomids);
-  molecule *mol = World::getInstance().
-      getMolecule(MoleculeById(atoms[0]->getMolecule()->getId()));
-  if (mol != atoms[1]->getMolecule()) {
-    STATUS("The two selected atoms must belong to the same molecule.");
-    return Action::failure;
-  }
-  // gather undo information
-  const double olddistance = atoms[0]->getPosition().distance(atoms[1]->getPosition());
-  const double newdistance = params.bonddistance.get();
-  LOG(1, "INFO: Old bond distance is " << olddistance << ", stretching to " << newdistance << ".");
+
+  const Vector NormalVector = (atoms[0]->getPosition() - atoms[1]->getPosition())* (1./olddistance);
+  const double shift = 0.5*(newdistance - olddistance);
+  std::vector<Vector> Shift(2);
+  Shift[0] = shift * NormalVector;
+  Shift[1] = -shift * NormalVector;
+  Box &domain = World::getInstance().getDomain();
 
   // Assume the selected bond splits the molecule into two parts, each one on
@@ -133 +114 @@
   // leaving the other positions untouched.
 
-  // convert BondGraph into boost::graph
+  // get nodes on either side of selected bond via BFS discovery
   BoostGraphCreator BGcreator;
   BGcreator.createFromMolecule(*mol,
       boost::bind(addEdgePredicate, _1, boost::ref(atomids)));
-  const BoostGraphCreator::UndirectedGraph &molgraph = BGcreator.get();
-
-  const size_t num_vertices = BGcreator.getNumVertices();
-  std::vector< std::vector<size_t> > distances;
-  for (size_t i=0;i<2;++i) {
-    distances.push_back(std::vector<size_t>(num_vertices, num_vertices+1)); // set distance to num+1
-    distances[i][atomids[i]] = 0;
-    boost::breadth_first_search(
-        molgraph,
-        boost::vertex(atomids[i], molgraph),
-        boost::visitor(record_distance(&(distances[i][0]))));
-    LOG(3, "DEBUG: From atom #" << atomids[i]
-        << " BFS discovered the following distances " << distances[i]);
-  }
-
-  const Vector NormalVector = (atoms[0]->getPosition() - atoms[1]->getPosition())* (1./olddistance);
-  const double shift = 0.5*(newdistance - olddistance);
-  std::vector<Vector> Shift(2);
-  Shift[0] = shift * NormalVector;
-  Shift[1] = -shift * NormalVector;
-  Box &domain = World::getInstance().getDomain();
-  std::vector< std::vector<size_t> > bondside_sets(2);
-
-  // Check whether there are common nodes in each set of distances
-  for (size_t i=0;i<num_vertices;++i)
-    if ((distances[0][i] != (num_vertices+1)) && (distances[1][i] != (num_vertices+1))) {
-      ELOG(2, "Node #" << i << " is reachable from either side of bond, hence must by cyclic."
+  BreadthFirstSearchGatherer NodeGatherer(BGcreator);
+  std::vector< BoostGraphHelpers::Nodeset_t > bondside_sets(2);
+  for(size_t j=0;j<2;++j) {
+    bondside_sets[j] = NodeGatherer(atoms[j]->getId());
+    std::sort(bondside_sets[j].begin(), bondside_sets[j].end());
+  }
+
+  // simple test whether bond has split the system in two disjoint sets or not
+  bool isCyclic = false;
+  if ((bondside_sets[0].size() + bondside_sets[1].size()) > BGcreator.getNumVertices()) {
+      // Check whether there are common nodes in each set of distances
+    if (BoostGraphHelpers::isCommonNodeInVector(bondside_sets[0], bondside_sets[1])) {
+      ELOG(2, "Sets contain common node, hence bond must have been by cyclic."
           << " Shifting only bond partners.");
       for(size_t j=0;j<2;++j) {
-        bondside_sets[j].push_back(atoms[j]->getId());
+        bondside_sets[j].clear();
+        bondside_sets[j].push_back( atomids[j] );
         const Vector &position = atoms[j]->getPosition();
         atoms[j]->setPosition( domain.enforceBoundaryConditions(position+Shift[j]) );
       }
-      break;
-    }
+      isCyclic = true;
+    }
+  }
 
   // go through the molecule and stretch each atom in either set of nodes
-  if (bondside_sets[0].empty()) {
-    const BoostGraphCreator::index_map_t index_map = BGcreator.getIndexMap();
+  if (!isCyclic) {
     for (molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
       const Vector &position = (*iter)->getPosition();
       // for each atom determine in which set of nodes it is and shift accordingly
-      size_t i=0;
-      const size_t nodeindex = boost::get(index_map, boost::vertex((*iter)->getId(), molgraph));
-      for (;i<2;++i) {
-        if (distances[i][nodeindex] != (num_vertices+1)) {
-          (*iter)->setPosition( domain.enforceBoundaryConditions(position+Shift[i]) );
-          bondside_sets[i].push_back((*iter)->getId());
-          break;
-        }
-      }
-      if (i==2) {
+      const atomId_t &atomid = (*iter)->getId();
+      if (std::binary_search(bondside_sets[0].begin(), bondside_sets[0].end(), atomid)) {
+        (*iter)->setPosition( domain.enforceBoundaryConditions(position+Shift[0]) );
+      } else if (std::binary_search(bondside_sets[1].begin(), bondside_sets[1].end(), atomid)) {
+        (*iter)->setPosition( domain.enforceBoundaryConditions(position+Shift[1]) );
+      } else {
         ELOG(1, "Atom " << *iter << " is not contained on either side of bond? Undoing done shifts");
         // Have to undo shifts
-        for (i=0;i<2;++i) {
-          for (std::vector<size_t>::const_iterator iter = bondside_sets[i].begin();
+        for (size_t i=0;i<2;++i) {
+          for (BoostGraphHelpers::Nodeset_t::const_iterator iter = bondside_sets[i].begin();
               iter != bondside_sets[i].end(); ++iter) {
             atom &walker = *World::getInstance().getAtom(AtomById(*iter));
@@ -214 +178 @@
   Box &domain = World::getInstance().getDomain();
   for (size_t i=0;i<2;++i) {
-    for (std::vector<size_t>::const_iterator iter = state->bondside_sets[i].begin();
+    for (BoostGraphHelpers::Nodeset_t::const_iterator iter = state->bondside_sets[i].begin();
         iter != state->bondside_sets[i].end(); ++iter) {
       atom &walker = *World::getInstance().getAtom(AtomById(*iter));
@@ -230 +194 @@
   Box &domain = World::getInstance().getDomain();
   for (size_t i=0;i<2;++i) {
-    for (std::vector<size_t>::const_iterator iter = state->bondside_sets[i].begin();
+    for (BoostGraphHelpers::Nodeset_t::const_iterator iter = state->bondside_sets[i].begin();
         iter != state->bondside_sets[i].end(); ++iter) {
       atom &walker = *World::getInstance().getAtom(AtomById(*iter));

src/Actions/MoleculeAction/StretchBondAction.def

@@ -21 +21 @@
 (RangeValidator< double >(0., 10.))
 
-#define statetypes (std::vector<Vector>)(std::vector< std::vector<size_t> >)(const molecule *)
+#define statetypes (std::vector<Vector>)(std::vector< BoostGraphHelpers::Nodeset_t >)(const molecule *)
 #define statereferences (Shift)(bondside_sets)(mol)
 

src/Graph/BoostGraphCreator.cpp

@@ -85 +85 @@
       std::back_inserter(atomids), getAtomId);
   ASSERT( _atoms.size() == atomids.size(),
-      "BoostGraphCreator::createFromAtom() - atomids and atoms differ in size?");
+      "BoostGraphCreator::createFromAtom() - atomids "
+      +toString(atomids.size())+" and atoms "+toString(_atoms.size())
+      +" differ in size?");
   std::sort(atomids.begin(), atomids.end());
   const predicate_t predicate = boost::bind(inSetPredicate, boost::ref(atomids), _1);
@@ -95 +97 @@
 }
 
+BoostGraphCreator::nodeId_t BoostGraphCreator::getNodeId(
+    const atomId_t &_atomid) const
+{
+  atomids_nodeids_t::const_iterator iter =
+      atomids_nodeids.find(_atomid);
+  return (iter == atomids_nodeids.end()) ? (nodeId_t)-1 : iter->second;
+}
+
 template <typename iterator>
 void BoostGraphCreator::createFromRange(
@@ -101 +111 @@
     const size_t &_no_nodes,
     const predicate_t &_pred
-    ) {
-  // convert BondGraph into boost::graph
-  UndirectedGraph molgraph(_no_nodes);
+    )
+{
+  graph = UndirectedGraph();
+
+  // add vertices
   for(iterator iter = _begin; iter != _end; ++iter) {
-    LOG(2, "DEBUG: Looking at node " << (*iter)->getId());
+    const atomId_t atomid = (*iter)->getId();
+    Vertex v = boost::add_vertex(atomid, graph);
+    const atomId_t vertexname = boost::get(boost::get(boost::vertex_name, graph), v);
+    const nodeId_t vertexindex = boost::get(boost::get(boost::vertex_index, graph), v);
+    LOG(2, "DEBUG: Adding node " << vertexindex << " associated to atom #" << vertexname);
+    ASSERT( vertexname == atomid,
+        "BoostGraphCreator::createFromRange() - atomid "+toString(atomid)
+        +" is not name of vertex "+toString(vertexname)+".");
+    atomids_nodeids.insert( std::make_pair(vertexname, vertexindex) );
+  }
+
+  // add edges
+  for(iterator iter = _begin; iter != _end; ++iter) {
+    LOG(2, "DEBUG: Looking at atom " << (*iter)->getId());
     const BondList& ListOfBonds = (*iter)->getListOfBonds();
     for(BondList::const_iterator bonditer = ListOfBonds.begin();
         bonditer != ListOfBonds.end(); ++bonditer) {
+      LOG(2, "DEBUG: Looking at bond " << *(*bonditer));
       const atomId_t leftid = (*bonditer)->leftatom->getId();
+      const nodeId_t leftnodeid = getNodeId(leftid);
       const atomId_t rightid = (*bonditer)->rightatom->getId();
+      const nodeId_t rightnodeid = getNodeId(rightid);
       // only pick each bond once and evaluate predicate
       if ((leftid == (*iter)->getId())
           && (_pred(**bonditer))) {
-        LOG(3, "DEBUG: ADDING edge " << leftid << " <-> " << rightid);
-        boost::add_edge(leftid, rightid, graph);
+        LOG(3, "DEBUG: ADDING edge " << leftnodeid << " <-> " << rightnodeid);
+        boost::add_edge(leftnodeid, rightnodeid, graph);
       } else {
-        LOG(3, "DEBUG: Discarding edge " << leftid << " <-> " << rightid);
+        LOG(3, "DEBUG: Discarding edge " << leftnodeid << " <-> " << rightnodeid);
       }
     }

src/Graph/BoostGraphCreator.hpp

@@ -15 +15 @@
 #endif
 
+#include <map>
 #include <vector>
 
 #include <boost/function.hpp>
 #include <boost/graph/adjacency_list.hpp>
+
+#include "types.hpp"
 
 class atom;
@@ -32 +35 @@
   //!> typedef for an undirected graph using boost::graph
   typedef boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS,
-      boost::no_property, boost::no_property > UndirectedGraph;
+      boost::property<boost::vertex_name_t, atomId_t>, boost::no_property > UndirectedGraph;
   //!> typedef for a map of graph node indices
-  typedef boost::property_map < boost::adjacency_list <>, boost::vertex_index_t >::type index_map_t;
+  typedef boost::property_map < UndirectedGraph, boost::vertex_index_t >::type index_map_t;
+  typedef boost::property_map < UndirectedGraph, boost::vertex_index_t >::const_type const_index_map_t;
+  //!> typedef for a map of graph node indices
+  typedef boost::property_map < UndirectedGraph, boost::vertex_name_t >::type name_map_t;
+  typedef boost::property_map < UndirectedGraph, boost::vertex_name_t >::const_type const_name_map_t;
   //!> typedef for the  predicate to evaluate for adding the current edge or not
   typedef boost::function<bool (const bond &)> predicate_t;
+  //!> typedef for a Vertex
+  typedef boost::graph_traits<UndirectedGraph>::vertex_descriptor Vertex;
+  //!> typedef for vertex iterator
+  typedef boost::graph_traits<UndirectedGraph>::vertex_iterator vertex_iter;
 
+  //!> typedef for a node id
+  typedef size_t nodeId_t;
+  //!> typedef for map converting between node id in graph and the associated atomic id
+  typedef std::map<atomId_t, nodeId_t> atomids_nodeids_t;
 
   /** Creates the boost::graph using all atoms and bonds in the given \a _mol.
@@ -89 +104 @@
   }
 
+  /** Returns the node id to a given atom id \a _atomid.
+   *
+   * \param _atomid atom id
+   * \return node id
+   */
+  nodeId_t getNodeId(const atomId_t &_atomid) const;
+
 private:
   /** General purpose function that contains the internal logic of walking the
@@ -114 +136 @@
   //!> internal graph that is created by creator functions
   UndirectedGraph graph;
+  //!> external property map for all the atomic ids of each graph node
+  atomids_nodeids_t atomids_nodeids;
 };
 

src/Graph/Makefile.am

@@ -3 +3 @@
 
 GRAPHSOURCE = \
+        Graph/AdjacencyList.cpp \
         Graph/BondGraph.cpp \
         Graph/BoostGraphCreator.cpp \
+        Graph/BreadthFirstSearchGatherer.cpp \
         Graph/BuildInducedSubgraph.cpp \
-        Graph/AdjacencyList.cpp \
         Graph/ConnectedSubgraph.cpp \
         Graph/CyclicStructureAnalysis.cpp \
@@ -12 +13 @@
                                   
 GRAPHHEADER = \
+        Graph/AdjacencyList.hpp \
         Graph/BondGraph.hpp \
         Graph/BoostGraphCreator.hpp \
+        Graph/BoostGraphHelpers.hpp \
+        Graph/BreadthFirstSearchGatherer.hpp \
         Graph/BuildInducedSubgraph.hpp \
-        Graph/AdjacencyList.hpp \
         Graph/ConnectedSubgraph.hpp \
         Graph/CyclicStructureAnalysis.hpp \