Context Navigation

-              rf5ea10
+              r99c705
 #include <iostream>
+#include <list>
+#include <map>
 #include <string>
 #include <vector>
+#include <boost/graph/copy.hpp>
 #include "CodePatterns/Log.hpp"
 …
 #include "Action_impl_pre.hpp"
 /** =========== define the function ====================== */
+template<class map_t, class key_t>
+static int getElementFromMap(
+    const map_t &_map,
+    const key_t &_keyname)
+{
+  typename map_t::const_iterator iter = _map.find(_keyname);
+  ASSERT( iter != _map.end(),
+      "getElementFromMap() -  cannot find key "+toString(_keyname)+" in map.");
+  return iter->second;
+}
+typedef std::pair<BoostGraphCreator::Vertex, BoostGraphCreator::Vertex> edge_by_vertices_t;
+typedef std::map<edge_by_vertices_t, int> edge_valence_t;
+typedef std::vector<edge_by_vertices_t> edges_t;
+typedef std::vector<int> degrees_t;
+typedef std::list< degrees_t > list_of_edge_degrees_t;
+typedef std::map< unsigned int, int> type_valence_map_t;
+typedef std::map< edge_by_vertices_t , size_t> edge_index_t;
+static int powerset_edge_degree(
+    list_of_edge_degrees_t &list_of_edge_degrees,
+    degrees_t &current_degrees,
+    const size_t remaining_edges,
+    const size_t current_index,
+    const edges_t &edges,
+    const edge_valence_t &edge_valence)
+{
+  if (remaining_edges == 0) {
+    list_of_edge_degrees.push_back(current_degrees);
+    LOG(3, "DEBUG: Adding [" << current_degrees << "] list of degrees.");
+    return 1;
+  } else {
+    size_t no_combinations = 0;
+    // get valence
+    const edge_by_vertices_t &e = edges[current_index];
+    const int current_max_degree = getElementFromMap<
+        edge_valence_t,
+        edge_by_vertices_t>(edge_valence, e);
+    // create all combinations and recurse
+    for(current_degrees[current_index] = 1;
+        current_degrees[current_index] <= current_max_degree;
+        ++current_degrees[current_index]) {
+      no_combinations += powerset_edge_degree(
+          list_of_edge_degrees, current_degrees,
+          remaining_edges-1, current_index+1,
+          edges, edge_valence);
+    }
+    return no_combinations;
+  }
+}
+static int getValenceForVertex(
+    const BoostGraphCreator::name_map_t &name_map,
+    const Extractors::type_index_lookup_t &type_index_lookup,
+    const type_valence_map_t &type_valence_map,
+    const BoostGraphCreator::Vertex &v
+    )
+{
+  const atomId_t atomid_v = boost::get(name_map, v);
+  const Extractors::type_index_lookup_t::left_const_iterator iter_v = type_index_lookup.left.find(atomid_v);
+  ASSERT( iter_v != type_index_lookup.left.end(),
+      "getValenceForVertex() - cannot find "+toString(atomid_v));
+  const unsigned int elementno_v = iter_v->second;
+  return getElementFromMap<type_valence_map_t, unsigned int>(type_valence_map, elementno_v);
+}
+static edge_by_vertices_t getEdgeByVertices(
+    const BoostGraphCreator::Vertex& u,
+    const BoostGraphCreator::Vertex& v)
+{
+  if (u < v)
+    return std::make_pair(u,v);
+  else
+    return std::make_pair(v,u);
+}
+static edge_by_vertices_t getEdgeVerticesByEdge(
+    const BoostGraphCreator::Edge& e,
+    const BoostGraphCreator::UndirectedGraph &graph)
+{
+  const BoostGraphCreator::Vertex u = source(e, graph);
+  const BoostGraphCreator::Vertex v = source(e, graph);
+  return getEdgeByVertices(u,v);
+}
+struct SaturatedGraph
+{
+  //!> states whether each node fulfills the octet rule, i.e. has bonds weight by their degree equal to its valence
+  bool ValencesAllFulfilled;
+  //!> graph saturated with additional hydrogens
+  BoostGraphCreator::UndirectedGraph saturated_graph;
+  //!> type per edge in the graph
+  Extractors::type_index_lookup_t type_index_lookup;
+};
+template <class graph_type>
+static SaturatedGraph saturateGraph(
+    const graph_type &graph,
+    const Extractors::nodes_t &nodes,
+    const BoostGraphCreator::name_map_t &name_map,
+    const Extractors::type_index_lookup_t &type_index_lookup,
+    const type_valence_map_t &type_valence_map,
+    const edge_index_t &edge_index,
+    const edges_t &edges,
+    const degrees_t &degrees
+)
+{
+  SaturatedGraph newgraph;
+  boost::copy_graph(graph, newgraph.saturated_graph);
+  newgraph.ValencesAllFulfilled = true;
+  // need to translate type_index_lookup as Extractors::createHomologyKey..() only
+  // looks at the vertex index, not its name. When extracting the subgraph, we
+  // might extract node #1, but insert it as new node index #0. Therefore, we have
+  // to translate all indices in the lookup
+  const Extractors::index_map_t index_map = boost::get(boost::vertex_index, graph);
+  {
+    BoostGraphCreator::vertex_iter vp, vpend;
+    for (boost::tie(vp, vpend) = vertices(graph); vp != vpend; ++vp) {
+      BoostGraphCreator::Vertex v = *vp;
+      const atomId_t vindex = boost::get(name_map, v);
+      Extractors::type_index_lookup_t::left_const_iterator iter = type_index_lookup.left.find(vindex);
+      if (iter != type_index_lookup.left.end()) {
+        const Extractors::node_t nodeid = boost::get(index_map, v);
+        newgraph.type_index_lookup.left.insert( std::make_pair(nodeid, iter->second) );
+        LOG(4, "DEBUG: Adding type to index lookup for vertex " << iter->first);
+      }
+    }
+  }
+  size_t no_nodes = boost::num_vertices(newgraph.saturated_graph);
+  // step through every node, sum up its degrees
+  BoostGraphCreator::vertex_iter vp, vpend;
+  for (boost::tie(vp, vpend) = vertices(graph); vp != vpend; ++vp) {
+    BoostGraphCreator::Vertex v = *vp;
+    const atomId_t vindex = boost::get(name_map, v);
+    LOG(3, "DEBUG: Current vertex is #" << vindex);
+    const int max_valence = getValenceForVertex(name_map, type_index_lookup, type_valence_map, v);
+    int total_valence = 0;
+    typename boost::graph_traits<graph_type>::out_edge_iterator i, end;
+    for (boost::tie(i, end) = boost::out_edges(v, graph); i != end; ++i) {
+      const BoostGraphCreator::Edge &e = *i;
+      const BoostGraphCreator::Vertex e1 = source(e, graph);
+      const BoostGraphCreator::Vertex e2 = target(e, graph);
+      const edge_by_vertices_t edge_by_vertices = getEdgeByVertices(e1,e2);
+      LOG(3, "DEBUG: Current edge is (" << edge_by_vertices << ")");
+      const size_t index = getElementFromMap<edge_index_t, edge_by_vertices_t>(
+          edge_index, edge_by_vertices);
+      ASSERT( index < edges.size(),
+          "saturateGraph() - index "+toString(index)+" out of range [0, "
+          +toString(edges.size())+"]");
+      total_valence += degrees[index];
+    }
+    LOG(3, "DEBUG: Vertex #" << vindex << " has total valence of "
+        << total_valence << " and " << max_valence << " max valence.");
+    // add hydrogens till valence is depleted
+    newgraph.ValencesAllFulfilled &= (total_valence <= max_valence);
+    for (int remaining_valence = total_valence; remaining_valence < max_valence;
+        ++remaining_valence) {
+      // add hydrogen node to this node
+      LOG(4, "DEBUG: Adding node " << no_nodes << " with type " << Extractors::ParticleType_t(1));
+      newgraph.type_index_lookup.left.insert(
+          std::make_pair(no_nodes, Extractors::ParticleType_t(1)));
+      BoostGraphCreator::Vertex h = boost::add_vertex(no_nodes, newgraph.saturated_graph);
+      ++no_nodes;
+      // add edge to hydrogen
+      std::pair<BoostGraphCreator::Edge, bool> edge_inserter =
+          boost::add_edge(v, h, newgraph.saturated_graph);
+    }
+    LOG(2, "DEBUG: Added " << (max_valence - total_valence)
+        << " hydrogens to vertex #" << vindex);
+  }
+  return newgraph;
+}
+template <typename graph_t>
+BoostGraphCreator::UndirectedGraph extractSubgraph(
+    const graph_t &_graph,
+    const Extractors::nodes_t &nodes)
+{
+  BoostGraphCreator::UndirectedGraph subgraph;
+  const Extractors::index_map_t index_map = boost::get(boost::vertex_index, _graph);
+  // add nodes
+  BoostGraphCreator::vertex_iter viter, vend;
+  for (boost::tie(viter, vend) = boost::vertices(_graph); viter != vend; ++viter) {
+    const size_t nodeid = boost::get(index_map, *viter);
+    if (nodes.find(nodeid) != nodes.end()) {
+      boost::add_vertex(*viter, subgraph);
+      LOG(4, "DEBUG: Adding node " << *viter << " to subgraph.");
+    }
+  }
+  // add edges
+  for (boost::tie(viter, vend) = boost::vertices(_graph); viter != vend; ++viter) {
+    BoostGraphCreator::UndirectedGraph::in_edge_iterator i, end;
+    for (boost::tie(i, end) = boost::in_edges(boost::vertex(*viter, _graph), _graph);
+        i != end; ++i) {
+      const Extractors::node_t sourceindex = boost::get(index_map, boost::source(*i, _graph));
+      const Extractors::node_t targetindex = boost::get(index_map, boost::target(*i, _graph));
+      if ((nodes.find(sourceindex) != nodes.end()) && (nodes.find(targetindex) != nodes.end())
+          && (sourceindex < targetindex)) {
+        boost::add_edge(sourceindex, targetindex, subgraph);
+        LOG(4, "DEBUG: Adding edge " << sourceindex << "<->" << targetindex << " to subgraph.");
+      }
+    }
+  }
+  return subgraph;
+}
 ActionState::ptr GraphChemicalSpaceEvaluatorAction::performCall() {
   // create boost::graph from graph6 string
+  /// 1. create boost::graph from graph6 string
   BoostGraphCreator BGC;
   BGC.createFromGraph6String(params.graph_string.get());
 …
   Extractors::index_map_t index_map = boost::get(boost::vertex_index, graph);
   LOG(2, "DEBUG: We have " << boost::num_vertices(graph) << " nodes in the fragment graph.");
   // associate type with a node index
+  LOG(1, "INFO: We have " << boost::num_vertices(graph) << " nodes in the fragment graph.");
+  /// 2. associate type with a node index
   Extractors::type_index_lookup_t type_index_lookup;
   std::vector<unsigned int> elementnumbers;
   elementnumbers.reserve(params.graph_elements.get().size());
+  size_t i = 0;
   for (std::vector<const element *>::const_iterator iter = params.graph_elements.get().begin();
       iter != params.graph_elements.get().end(); ++iter) {
     elementnumbers.push_back((*iter)->getAtomicNumber());
+    LOG(3, "DEBUG: Adding [" << i << ", " << elementnumbers.back() << "] type to index lookup.");
     type_index_lookup.left.insert(
+        std::make_pair(
+            std::distance(iter, params.graph_elements.get().begin()),
+            Extractors::ParticleType_t(elementnumbers.back())));
+  }
+  // create all induced subgraphs
+        std::make_pair(i++, Extractors::ParticleType_t(elementnumbers.back())));
+  }
+  if (i != boost::num_vertices(graph)) {
+    ELOG(1, "Not enough elements given.");
+    return Action::failure;
+  }
+  /// 3a. get range of valence for each node
+  type_valence_map_t type_valence_map;
+  for (std::vector<const element *>::const_iterator iter = params.graph_elements.get().begin();
+      iter != params.graph_elements.get().end(); ++iter) {
+    if (type_valence_map.count((*iter)->getAtomicNumber()) == 0) {
+      LOG(3, "DEBUG: Adding [" << (*iter)->getAtomicNumber()
+          << ", " << (*iter)->getNoValenceOrbitals() << "] type to valence lookup.");
+      type_valence_map.insert( std::make_pair(
+          (*iter)->getAtomicNumber(),
+          (*iter)->getNoValenceOrbitals()) );
+    }
+  }
+  // get map for retrieving vertex index for a given vertex
+  const BoostGraphCreator::name_map_t name_map = boost::get(boost::vertex_name, graph);
+  /// 3b. generate range of degree for each edge
+  edge_valence_t edge_valence;
+  edges_t edges;
+  edge_index_t edge_index;
+  edges.reserve(boost::num_edges(graph));
+  {
+    size_t index = 0;
+    BoostGraphCreator::edge_iter i, end;
+    for (boost::tie(i, end) = boost::edges(graph); i != end; ++i) {
+      static int MAX_DEGREE = 3;
+      const BoostGraphCreator::Edge &e = *i;
+      /// 4a. place all edges into a vector
+      const BoostGraphCreator::Vertex u = source(e, graph);
+      const BoostGraphCreator::Vertex v = target(e, graph);
+      const edge_by_vertices_t edge_by_vertices = getEdgeByVertices(u,v);
+      edges.push_back(edge_by_vertices);
+      edge_index.insert( std::make_pair(
+          edge_by_vertices, index++) );
+      const int valence_u = getValenceForVertex(name_map, type_index_lookup, type_valence_map, u);
+      const int valence_v = getValenceForVertex(name_map, type_index_lookup, type_valence_map, v);
+      const int max_degree = std::max(1,
+          std::min(
+              std::min(valence_u, MAX_DEGREE),
+              std::min(valence_v, MAX_DEGREE)
+              ));
+      edge_valence.insert(
+          std::make_pair(edge_by_vertices, max_degree)
+          );
+      LOG(3, "DEBUG: Adding [" << std::make_pair(u,v) << ", "
+          << max_degree << "] edge to valence lookup.");
+    }
+  }
+  /// 4. for each graph with type per node and degree per edge set, add hydrogen nodes and edges
+  // go through powerset over all edges
+  /// 4b. have a recursive function that puts together the powerset over all max_degree
+  list_of_edge_degrees_t list_of_edge_degrees;
+  size_t no_combinations = 0;
+  {
+    size_t remaining_edges = edges.size();
+    size_t current_index = 0;
+    degrees_t current_degrees(edges.size(), 0);
+    no_combinations = powerset_edge_degree(
+        list_of_edge_degrees, current_degrees,
+        remaining_edges, current_index,
+        edges, edge_valence);
+  }
+  LOG(1, "INFO: Added " << no_combinations << " graph degree combinations.");
+  /// 5. dissect graph into a list of subgraphs
+  typedef std::vector<BoostGraphCreator::UndirectedGraph> subgraphs_t;
+  subgraphs_t allsubgraphs;
   const size_t nodes_in_graph = boost::num_vertices(graph);
   LOG(2, "DEBUG: There are " << nodes_in_graph << " nodes in the graph.");
 …
     LOG(3, "DEBUG: Nodes per level are " << nodes_per_level);
     /// 6. construct all possible induced connected subgraphs with this map (see fragmentation)
+    /// 5c. construct all possible induced connected subgraphs with this map (see fragmentation)
     Extractors::nodes_t nodes;
     // rootindex is always contained
 …
     Extractors::level_t level = 0;
+    Extractors::generateAllInducedConnectedSubgraphs(
+        nodes_in_graph-1, level, nodes, set_of_nodes, nodes_per_level, graph, distances, index_map);
+    for (size_t i=0;i<nodes_in_graph; ++i)
+      Extractors::generateAllInducedConnectedSubgraphs(
+          i, level, nodes, set_of_nodes, nodes_per_level, graph, distances, index_map);
     LOG(2, "DEBUG: We have found " << set_of_nodes.size() << " unique induced, connected subgraphs.");
+  }
+  typedef std::vector<
+      std::pair<
+          BoostGraphCreator::UndirectedGraph,
+          degrees_t> > graphs_with_degrees_t;
+  // 6. for every combination saturate fragments for lookup and energy contribution summation
   const HomologyContainer &container = World::getInstance().getHomologies();
+  /// go through each induced, connected subgraph
+  for (Extractors::set_of_nodes_t::const_iterator nodesiter = set_of_nodes.begin();
+      nodesiter != set_of_nodes.end(); ++nodesiter) {
+    /// convert its nodes into a HomologyGraph using the graph and elements (types for index)
+    const Extractors::nodes_t &nodes = *nodesiter;
+    const HomologyGraph nodes_graph =
+        Extractors::createHomologyGraphFromNodes(nodes, type_index_lookup, graph, index_map);
+    /// Query HomologyContainer for this HomologyGraph.
+    HomologyContainer::range_t range = container.getHomologousGraphs(nodes_graph);
+    if (range.first == range.second) {
+      // range is empty, the fragment is unknown
+      ELOG(1, "Cannot find fragment graph " << nodes_graph << " to graph " << elementnumbers);
+    } else {
+      // list first energy
+      LOG(1, "Fragment graph " << nodes_graph << " has energy " << range.first->second.energy);
+    }
+  for (list_of_edge_degrees_t::const_iterator degrees_iter = list_of_edge_degrees.begin();
+      degrees_iter != list_of_edge_degrees.end(); ++degrees_iter) {
+    const degrees_t &current_degrees = *degrees_iter;
+    // simply sum up the contributions from all fragments for the total energy
+    // of this particular combination, given the graph, the element for each node,
+    // and the bond degrees for each edge
+    double total_energy = 0.;
+    // get fragment
+    for (Extractors::set_of_nodes_t::const_iterator nodesiter = set_of_nodes.begin();
+        nodesiter != set_of_nodes.end(); ++nodesiter) {
+      const Extractors::nodes_t &current_nodes = *nodesiter;
+      LOG(2, "DEBUG: Nodes for graph are " << current_nodes);
+      // create subgraph
+      Extractors::UndirectedGraph newgraph = extractSubgraph(graph, current_nodes);
+      const BoostGraphCreator::name_map_t new_name_map =
+          boost::get(boost::vertex_name, newgraph);
+      {
+        BoostGraphCreator::vertex_iter viter, vend;
+        boost::tie(viter, vend) = boost::vertices(newgraph);
+        for (; viter != vend; ++viter) {
+          const atomId_t vindex = boost::get(new_name_map, *viter);
+          LOG(4, "DEBUG: Copied node has index " << vindex);
+        }
+      }
+      // saturate subgraph: add to nodes, add to type_index_lookup, index_map
+        // saturate any leftover valence with hydrogen nodes
+      SaturatedGraph fullgraph =
+          saturateGraph(newgraph, current_nodes, new_name_map,
+              type_index_lookup, type_valence_map, edge_index,
+              edges, current_degrees);
+      // don't add if one of the node's valence cannot be fulfilled
+      if (!fullgraph.ValencesAllFulfilled)
+        continue;
+      const Extractors::index_map_t subgraph_index_map =
+          boost::get(boost::vertex_index, fullgraph.saturated_graph);
+      const BoostGraphCreator::name_map_t subgraph_name_map =
+          boost::get(boost::vertex_name, fullgraph.saturated_graph);
+      Extractors::nodes_t saturated_nodes;
+      BoostGraphCreator::vertex_iter viter, vend;
+      boost::tie(viter, vend) = boost::vertices(fullgraph.saturated_graph);
+      for (; viter != vend; ++viter) {
+        const Extractors::node_t vindex = boost::get(subgraph_index_map, *viter);
+        LOG(4, "DEBUG: Adding node index " << vindex << " to saturated nodes list.");
+        saturated_nodes.insert(vindex);
+      }
+      /// convert its nodes into a HomologyGraph using the graph and elements (types for index)
+      const HomologyGraph nodes_graph =
+          Extractors::createHomologyGraphFromNodes(
+              saturated_nodes,
+              fullgraph.type_index_lookup,
+              fullgraph.saturated_graph,
+              subgraph_index_map);
+      LOG(2, "DEBUG: Looking for graph " << nodes_graph);
+      /// Query HomologyContainer for this HomologyGraph.
+      HomologyContainer::range_t range = container.getHomologousGraphs(nodes_graph);
+      if (range.first == range.second) {
+        // range is empty, the fragment is unknown
+        ELOG(1, "Cannot find fragment graph " << nodes_graph << " to graph " << elementnumbers);
+      } else {
+        // list first energy
+        LOG(1, "Fragment graph " << nodes_graph << " has contribution " << range.first->second.energy);
+        total_energy += range.first->second.energy;
+      }
+    }
+    LOG(1, "The graph with degrees " << *degrees_iter << " has a total energy of " << total_energy);
+  }

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 99c705 for src/Actions

Legend:

src/Actions/GraphAction/ChemicalSpaceEvaluatorAction.cpp

Download in other formats: