| 1 | /*
 | 
|---|
| 2 |  * Project: MoleCuilder
 | 
|---|
| 3 |  * Description: creates and alters molecular systems
 | 
|---|
| 4 |  * Copyright (C)  2017 Frederik Heber. All rights reserved.
 | 
|---|
| 5 |  * 
 | 
|---|
| 6 |  *
 | 
|---|
| 7 |  *   This file is part of MoleCuilder.
 | 
|---|
| 8 |  *
 | 
|---|
| 9 |  *    MoleCuilder is free software: you can redistribute it and/or modify
 | 
|---|
| 10 |  *    it under the terms of the GNU General Public License as published by
 | 
|---|
| 11 |  *    the Free Software Foundation, either version 2 of the License, or
 | 
|---|
| 12 |  *    (at your option) any later version.
 | 
|---|
| 13 |  *
 | 
|---|
| 14 |  *    MoleCuilder is distributed in the hope that it will be useful,
 | 
|---|
| 15 |  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 | 
|---|
| 16 |  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 | 
|---|
| 17 |  *    GNU General Public License for more details.
 | 
|---|
| 18 |  *
 | 
|---|
| 19 |  *    You should have received a copy of the GNU General Public License
 | 
|---|
| 20 |  *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 | 
|---|
| 21 |  */
 | 
|---|
| 22 | 
 | 
|---|
| 23 | /*
 | 
|---|
| 24 |  * ChemicalSpaceEvaluator.cpp
 | 
|---|
| 25 |  *
 | 
|---|
| 26 |  *  Created on: Sep 26, 2017
 | 
|---|
| 27 |  *      Author: heber
 | 
|---|
| 28 |  */
 | 
|---|
| 29 | 
 | 
|---|
| 30 | // include config.h
 | 
|---|
| 31 | #ifdef HAVE_CONFIG_H
 | 
|---|
| 32 | #include <config.h>
 | 
|---|
| 33 | #endif
 | 
|---|
| 34 | 
 | 
|---|
| 35 | //#include "CodePatterns/MemDebug.hpp"
 | 
|---|
| 36 | 
 | 
|---|
| 37 | #include <iostream>
 | 
|---|
| 38 | #include <list>
 | 
|---|
| 39 | #include <map>
 | 
|---|
| 40 | #include <string>
 | 
|---|
| 41 | #include <vector>
 | 
|---|
| 42 | 
 | 
|---|
| 43 | #include <boost/graph/copy.hpp>
 | 
|---|
| 44 | 
 | 
|---|
| 45 | #include "CodePatterns/Log.hpp"
 | 
|---|
| 46 | 
 | 
|---|
| 47 | #include "Actions/GraphAction/ChemicalSpaceEvaluatorAction.hpp"
 | 
|---|
| 48 | 
 | 
|---|
| 49 | #include "Element/element.hpp"
 | 
|---|
| 50 | #include "Graph/BoostGraphCreator.hpp"
 | 
|---|
| 51 | #include "Fragmentation/Homology/HomologyContainer.hpp"
 | 
|---|
| 52 | #include "FunctionApproximation/Extractors.hpp"
 | 
|---|
| 53 | #include "World.hpp"
 | 
|---|
| 54 | 
 | 
|---|
| 55 | #include <boost/graph/adjacency_list.hpp>
 | 
|---|
| 56 | 
 | 
|---|
| 57 | using namespace MoleCuilder;
 | 
|---|
| 58 | 
 | 
|---|
| 59 | // and construct the stuff
 | 
|---|
| 60 | #include "ChemicalSpaceEvaluatorAction.def"
 | 
|---|
| 61 | #include "Action_impl_pre.hpp"
 | 
|---|
| 62 | /** =========== define the function ====================== */
 | 
|---|
| 63 | template<class map_t, class key_t>
 | 
|---|
| 64 | static int getElementFromMap(
 | 
|---|
| 65 |     const map_t &_map,
 | 
|---|
| 66 |     const key_t &_keyname)
 | 
|---|
| 67 | {
 | 
|---|
| 68 |   typename map_t::const_iterator iter = _map.find(_keyname);
 | 
|---|
| 69 |   ASSERT( iter != _map.end(),
 | 
|---|
| 70 |       "getElementFromMap() -  cannot find key "+toString(_keyname)+" in map.");
 | 
|---|
| 71 |   return iter->second;
 | 
|---|
| 72 | }
 | 
|---|
| 73 | 
 | 
|---|
| 74 | typedef std::pair<BoostGraphCreator::Vertex, BoostGraphCreator::Vertex> edge_by_vertices_t;
 | 
|---|
| 75 | typedef std::map<edge_by_vertices_t, int> edge_valence_t;
 | 
|---|
| 76 | typedef std::vector<edge_by_vertices_t> edges_t;
 | 
|---|
| 77 | typedef std::vector<int> degrees_t;
 | 
|---|
| 78 | typedef std::list< degrees_t > list_of_edge_degrees_t;
 | 
|---|
| 79 | typedef std::map< unsigned int, int> type_valence_map_t;
 | 
|---|
| 80 | typedef std::map< edge_by_vertices_t , size_t> edge_index_t;
 | 
|---|
| 81 | 
 | 
|---|
| 82 | static int powerset_edge_degree(
 | 
|---|
| 83 |     list_of_edge_degrees_t &list_of_edge_degrees,
 | 
|---|
| 84 |     degrees_t ¤t_degrees,
 | 
|---|
| 85 |     const size_t remaining_edges,
 | 
|---|
| 86 |     const size_t current_index,
 | 
|---|
| 87 |     const edges_t &edges,
 | 
|---|
| 88 |     const edge_valence_t &edge_valence)
 | 
|---|
| 89 | {
 | 
|---|
| 90 |   if (remaining_edges == 0) {
 | 
|---|
| 91 |     list_of_edge_degrees.push_back(current_degrees);
 | 
|---|
| 92 |     LOG(3, "DEBUG: Adding [" << current_degrees << "] list of degrees.");
 | 
|---|
| 93 |     return 1;
 | 
|---|
| 94 |   } else {
 | 
|---|
| 95 |     size_t no_combinations = 0;
 | 
|---|
| 96 |     // get valence
 | 
|---|
| 97 |     const edge_by_vertices_t &e = edges[current_index];
 | 
|---|
| 98 |     const int current_max_degree = getElementFromMap<
 | 
|---|
| 99 |         edge_valence_t,
 | 
|---|
| 100 |         edge_by_vertices_t>(edge_valence, e);
 | 
|---|
| 101 |     // create all combinations and recurse
 | 
|---|
| 102 |     for(current_degrees[current_index] = 1;
 | 
|---|
| 103 |         current_degrees[current_index] <= current_max_degree;
 | 
|---|
| 104 |         ++current_degrees[current_index]) {
 | 
|---|
| 105 |       no_combinations += powerset_edge_degree(
 | 
|---|
| 106 |           list_of_edge_degrees, current_degrees,
 | 
|---|
| 107 |           remaining_edges-1, current_index+1,
 | 
|---|
| 108 |           edges, edge_valence);
 | 
|---|
| 109 |     }
 | 
|---|
| 110 |     return no_combinations;
 | 
|---|
| 111 |   }
 | 
|---|
| 112 | }
 | 
|---|
| 113 | 
 | 
|---|
| 114 | static int getValenceForVertex(
 | 
|---|
| 115 |     const BoostGraphCreator::name_map_t &name_map,
 | 
|---|
| 116 |     const Extractors::type_index_lookup_t &type_index_lookup,
 | 
|---|
| 117 |     const type_valence_map_t &type_valence_map,
 | 
|---|
| 118 |     const BoostGraphCreator::Vertex &v
 | 
|---|
| 119 |     )
 | 
|---|
| 120 | {
 | 
|---|
| 121 |   const atomId_t atomid_v = boost::get(name_map, v);
 | 
|---|
| 122 |   const Extractors::type_index_lookup_t::left_const_iterator iter_v = type_index_lookup.left.find(atomid_v);
 | 
|---|
| 123 |   ASSERT( iter_v != type_index_lookup.left.end(),
 | 
|---|
| 124 |       "getValenceForVertex() - cannot find "+toString(atomid_v));
 | 
|---|
| 125 |   const unsigned int elementno_v = iter_v->second;
 | 
|---|
| 126 |   return getElementFromMap<type_valence_map_t, unsigned int>(type_valence_map, elementno_v);
 | 
|---|
| 127 | }
 | 
|---|
| 128 | 
 | 
|---|
| 129 | static edge_by_vertices_t getEdgeByVertices(
 | 
|---|
| 130 |     const BoostGraphCreator::Vertex& u,
 | 
|---|
| 131 |     const BoostGraphCreator::Vertex& v)
 | 
|---|
| 132 | {
 | 
|---|
| 133 |   if (u < v)
 | 
|---|
| 134 |     return std::make_pair(u,v);
 | 
|---|
| 135 |   else
 | 
|---|
| 136 |     return std::make_pair(v,u);
 | 
|---|
| 137 | }
 | 
|---|
| 138 | 
 | 
|---|
| 139 | static edge_by_vertices_t getEdgeVerticesByEdge(
 | 
|---|
| 140 |     const BoostGraphCreator::Edge& e,
 | 
|---|
| 141 |     const BoostGraphCreator::UndirectedGraph &graph)
 | 
|---|
| 142 | {
 | 
|---|
| 143 |   const BoostGraphCreator::Vertex u = source(e, graph);
 | 
|---|
| 144 |   const BoostGraphCreator::Vertex v = target(e, graph);
 | 
|---|
| 145 |   return getEdgeByVertices(u,v);
 | 
|---|
| 146 | }
 | 
|---|
| 147 | 
 | 
|---|
| 148 | 
 | 
|---|
| 149 | struct SaturatedGraph
 | 
|---|
| 150 | {
 | 
|---|
| 151 |   //!> states whether each node fulfills the octet rule, i.e. has bonds weight by their degree equal to its valence
 | 
|---|
| 152 |   bool ValencesAllFulfilled;
 | 
|---|
| 153 |   //!> graph saturated with additional hydrogens
 | 
|---|
| 154 |   BoostGraphCreator::UndirectedGraph saturated_graph;
 | 
|---|
| 155 |   //!> type per edge in the graph
 | 
|---|
| 156 |   Extractors::type_index_lookup_t type_index_lookup;
 | 
|---|
| 157 | };
 | 
|---|
| 158 | 
 | 
|---|
| 159 | template <class graph_type>
 | 
|---|
| 160 | static SaturatedGraph saturateGraph(
 | 
|---|
| 161 |     const graph_type &graph,
 | 
|---|
| 162 |     const BoostGraphCreator::name_map_t &name_map,
 | 
|---|
| 163 |     const Extractors::type_index_lookup_t &type_index_lookup,
 | 
|---|
| 164 |     const type_valence_map_t &type_valence_map,
 | 
|---|
| 165 |     const edge_index_t &edge_index,
 | 
|---|
| 166 |     const edges_t &edges,
 | 
|---|
| 167 |     const degrees_t °rees
 | 
|---|
| 168 | )
 | 
|---|
| 169 | {
 | 
|---|
| 170 |   SaturatedGraph newgraph;
 | 
|---|
| 171 |   boost::copy_graph(graph, newgraph.saturated_graph);
 | 
|---|
| 172 |   newgraph.ValencesAllFulfilled = true;
 | 
|---|
| 173 | 
 | 
|---|
| 174 |   // need to translate type_index_lookup as Extractors::createHomologyKey..() only
 | 
|---|
| 175 |   // looks at the vertex index, not its name. When extracting the subgraph, we
 | 
|---|
| 176 |   // might extract node #1, but insert it as new node index #0. Therefore, we have
 | 
|---|
| 177 |   // to translate all indices in the lookup
 | 
|---|
| 178 |   const Extractors::index_map_t index_map = boost::get(boost::vertex_index, graph);
 | 
|---|
| 179 |   {
 | 
|---|
| 180 |     BoostGraphCreator::vertex_iter vp, vpend;
 | 
|---|
| 181 |     for (boost::tie(vp, vpend) = vertices(graph); vp != vpend; ++vp) {
 | 
|---|
| 182 |       BoostGraphCreator::Vertex v = *vp;
 | 
|---|
| 183 |       const atomId_t vindex = boost::get(name_map, v);
 | 
|---|
| 184 |       Extractors::type_index_lookup_t::left_const_iterator iter = type_index_lookup.left.find(vindex);
 | 
|---|
| 185 |       if (iter != type_index_lookup.left.end()) {
 | 
|---|
| 186 |         const Extractors::node_t nodeid = boost::get(index_map, v);
 | 
|---|
| 187 |         newgraph.type_index_lookup.left.insert( std::make_pair(nodeid, iter->second) );
 | 
|---|
| 188 |         LOG(4, "DEBUG: Adding type to index lookup for vertex " << iter->first);
 | 
|---|
| 189 |       }
 | 
|---|
| 190 |     }
 | 
|---|
| 191 |   }
 | 
|---|
| 192 | 
 | 
|---|
| 193 |   LOG(2, "DEBUG: Saturating graph with " << boost::num_vertices(newgraph.saturated_graph) << " nodes.");
 | 
|---|
| 194 |   size_t no_nodes = boost::num_vertices(newgraph.saturated_graph);
 | 
|---|
| 195 |   // step through every node, sum up its degrees
 | 
|---|
| 196 |   BoostGraphCreator::vertex_iter vp, vpend;
 | 
|---|
| 197 |   for (boost::tie(vp, vpend) = vertices(graph); vp != vpend; ++vp) {
 | 
|---|
| 198 |     BoostGraphCreator::Vertex v = *vp;
 | 
|---|
| 199 |     const atomId_t vindex = boost::get(name_map, v);
 | 
|---|
| 200 |     LOG(3, "DEBUG: Current vertex is #" << vindex);
 | 
|---|
| 201 |     const int max_valence = getValenceForVertex(name_map, type_index_lookup, type_valence_map, v);
 | 
|---|
| 202 |     int total_valence = 0;
 | 
|---|
| 203 |     typename boost::graph_traits<graph_type>::out_edge_iterator i, end;
 | 
|---|
| 204 |     for (boost::tie(i, end) = boost::out_edges(v, graph); i != end; ++i) {
 | 
|---|
| 205 |       const BoostGraphCreator::Edge &e = *i;
 | 
|---|
| 206 |       const edge_by_vertices_t edge_by_vertices = getEdgeVerticesByEdge(e, graph);
 | 
|---|
| 207 |       LOG(3, "DEBUG: Current edge is (" << edge_by_vertices << ")");
 | 
|---|
| 208 |       const size_t index = getElementFromMap<edge_index_t, edge_by_vertices_t>(
 | 
|---|
| 209 |           edge_index, edge_by_vertices);
 | 
|---|
| 210 |       ASSERT( index < edges.size(),
 | 
|---|
| 211 |           "saturateGraph() - index "+toString(index)+" out of range [0, "
 | 
|---|
| 212 |           +toString(edges.size())+"]");
 | 
|---|
| 213 |       total_valence += degrees[index];
 | 
|---|
| 214 |     }
 | 
|---|
| 215 |     LOG(3, "DEBUG: Vertex #" << vindex << " has total valence of "
 | 
|---|
| 216 |         << total_valence << " and " << max_valence << " max valence.");
 | 
|---|
| 217 |     // add hydrogens till valence is depleted
 | 
|---|
| 218 |     newgraph.ValencesAllFulfilled &= (total_valence <= max_valence);
 | 
|---|
| 219 |     for (int remaining_valence = total_valence; remaining_valence < max_valence;
 | 
|---|
| 220 |         ++remaining_valence) {
 | 
|---|
| 221 |       // add hydrogen node to this node
 | 
|---|
| 222 |       LOG(4, "DEBUG: Adding node " << no_nodes << " with type " << Extractors::ParticleType_t(1));
 | 
|---|
| 223 |       newgraph.type_index_lookup.left.insert(
 | 
|---|
| 224 |           std::make_pair(no_nodes, Extractors::ParticleType_t(1)));
 | 
|---|
| 225 |       BoostGraphCreator::Vertex h = boost::add_vertex(no_nodes, newgraph.saturated_graph);
 | 
|---|
| 226 |       ++no_nodes;
 | 
|---|
| 227 | 
 | 
|---|
| 228 |       // add edge to hydrogen
 | 
|---|
| 229 |       std::pair<BoostGraphCreator::Edge, bool> edge_inserter =
 | 
|---|
| 230 |           boost::add_edge(v, h, newgraph.saturated_graph);
 | 
|---|
| 231 |       ASSERT( edge_inserter.second,
 | 
|---|
| 232 |           "Failed to insert hydrogen into saturation graph.");
 | 
|---|
| 233 |     }
 | 
|---|
| 234 |     LOG(2, "DEBUG: Added " << (max_valence - total_valence)
 | 
|---|
| 235 |         << " hydrogens to vertex #" << vindex);
 | 
|---|
| 236 |   }
 | 
|---|
| 237 | 
 | 
|---|
| 238 |   return newgraph;
 | 
|---|
| 239 | }
 | 
|---|
| 240 | 
 | 
|---|
| 241 | template <typename graph_t>
 | 
|---|
| 242 | BoostGraphCreator::UndirectedGraph extractSubgraph(
 | 
|---|
| 243 |     const graph_t &_graph,
 | 
|---|
| 244 |     const Extractors::nodes_t &nodes,
 | 
|---|
| 245 |     const edge_index_t &_edge_index,
 | 
|---|
| 246 |     edge_index_t &_subgraph_edge_index)
 | 
|---|
| 247 | {
 | 
|---|
| 248 |   BoostGraphCreator::UndirectedGraph subgraph;
 | 
|---|
| 249 |   const Extractors::index_map_t index_map = boost::get(boost::vertex_index, _graph);
 | 
|---|
| 250 |   typedef std::map<Extractors::node_t, BoostGraphCreator::Vertex> graph_index_to_subgraph_vertex_t;
 | 
|---|
| 251 |   graph_index_to_subgraph_vertex_t graph_index_to_subgraph_vertex;
 | 
|---|
| 252 | 
 | 
|---|
| 253 |   // add nodes
 | 
|---|
| 254 |   BoostGraphCreator::vertex_iter viter, vend;
 | 
|---|
| 255 |   for (boost::tie(viter, vend) = boost::vertices(_graph); viter != vend; ++viter) {
 | 
|---|
| 256 |     const Extractors::node_t nodeid = boost::get(index_map, *viter);
 | 
|---|
| 257 |     if (nodes.find(nodeid) != nodes.end()) {
 | 
|---|
| 258 |       const BoostGraphCreator::Vertex vnew = boost::add_vertex(*viter, subgraph);
 | 
|---|
| 259 |       graph_index_to_subgraph_vertex.insert( std::make_pair(nodeid, vnew) );
 | 
|---|
| 260 |       LOG(4, "DEBUG: Adding node " << *viter << " to subgraph.");
 | 
|---|
| 261 |     }
 | 
|---|
| 262 |   }
 | 
|---|
| 263 |   const Extractors::index_map_t subgraph_index_map = boost::get(boost::vertex_index, subgraph);
 | 
|---|
| 264 | 
 | 
|---|
| 265 |   // add edges
 | 
|---|
| 266 |   for (boost::tie(viter, vend) = boost::vertices(_graph); viter != vend; ++viter) {
 | 
|---|
| 267 |     BoostGraphCreator::UndirectedGraph::in_edge_iterator i, end;
 | 
|---|
| 268 |     for (boost::tie(i, end) = boost::in_edges(boost::vertex(*viter, _graph), _graph);
 | 
|---|
| 269 |         i != end; ++i) {
 | 
|---|
| 270 |       const Extractors::node_t sourceindex = boost::get(index_map, boost::source(*i, _graph));
 | 
|---|
| 271 |       const Extractors::node_t targetindex = boost::get(index_map, boost::target(*i, _graph));
 | 
|---|
| 272 |       // we need to translate the vertex index from graph to subgraph
 | 
|---|
| 273 |       const Extractors::node_t subsourceindex = boost::get(
 | 
|---|
| 274 |           subgraph_index_map, graph_index_to_subgraph_vertex[sourceindex]);
 | 
|---|
| 275 |       const Extractors::node_t subtargetindex = boost::get(
 | 
|---|
| 276 |           subgraph_index_map, graph_index_to_subgraph_vertex[targetindex]);
 | 
|---|
| 277 |       if ((nodes.find(sourceindex) != nodes.end()) && (nodes.find(targetindex) != nodes.end())
 | 
|---|
| 278 |           && (sourceindex < targetindex)) {
 | 
|---|
| 279 |         // and we need to translate the edge index from the graph to the subgraph
 | 
|---|
| 280 |         const std::pair<BoostGraphCreator::Edge, bool> newedgeinserter =
 | 
|---|
| 281 |             boost::add_edge(subsourceindex, subtargetindex, subgraph);
 | 
|---|
| 282 |         ASSERT( newedgeinserter.second,
 | 
|---|
| 283 |             "extractSubgraph() - could not insert edge "+toString(subsourceindex)+"<->"
 | 
|---|
| 284 |             +toString(subtargetindex)+".");
 | 
|---|
| 285 |         const edge_by_vertices_t edge_by_vertices = getEdgeVerticesByEdge(*i, _graph);
 | 
|---|
| 286 |         const edge_index_t::const_iterator edgeiter = _edge_index.find(edge_by_vertices);
 | 
|---|
| 287 |         ASSERT( edgeiter != _edge_index.end(),
 | 
|---|
| 288 |             "extractSubgraph() - could not find edge "+toString(edge_by_vertices)+" in edge_index map." );
 | 
|---|
| 289 |         const edge_by_vertices_t subgraph_edge_by_vertices =
 | 
|---|
| 290 |             getEdgeVerticesByEdge(newedgeinserter.first, subgraph);
 | 
|---|
| 291 |         _subgraph_edge_index.insert( std::make_pair(subgraph_edge_by_vertices, edgeiter->second) );
 | 
|---|
| 292 |         LOG(4, "DEBUG: Adding edge " << sourceindex << "<->" << targetindex
 | 
|---|
| 293 |             << " in graph to subgraph as edge " << subsourceindex << "<->" << subtargetindex << ".");
 | 
|---|
| 294 |       }
 | 
|---|
| 295 |     }
 | 
|---|
| 296 |   }
 | 
|---|
| 297 | 
 | 
|---|
| 298 |   return subgraph;
 | 
|---|
| 299 | }
 | 
|---|
| 300 | 
 | 
|---|
| 301 | 
 | 
|---|
| 302 | bool compareEnergyContribution(
 | 
|---|
| 303 |       const std::pair<const HomologyGraph, HomologyContainer::value_t> &a,
 | 
|---|
| 304 |       const std::pair<const HomologyGraph, HomologyContainer::value_t> &b) {
 | 
|---|
| 305 |   return a.second.energy < b.second.energy;
 | 
|---|
| 306 | }
 | 
|---|
| 307 | 
 | 
|---|
| 308 | ActionState::ptr GraphChemicalSpaceEvaluatorAction::performCall() {
 | 
|---|
| 309 |   /// 1. create boost::graph from graph6 string
 | 
|---|
| 310 |   BoostGraphCreator BGC;
 | 
|---|
| 311 |   BGC.createFromGraph6String(params.graph_string.get());
 | 
|---|
| 312 | 
 | 
|---|
| 313 |   BoostGraphCreator::UndirectedGraph graph = BGC.get();
 | 
|---|
| 314 | 
 | 
|---|
| 315 |   Extractors::index_map_t index_map = boost::get(boost::vertex_index, graph);
 | 
|---|
| 316 |   LOG(1, "INFO: We have " << boost::num_vertices(graph) << " nodes in the fragment graph.");
 | 
|---|
| 317 | 
 | 
|---|
| 318 |   /// 2. associate type with a node index
 | 
|---|
| 319 |   Extractors::type_index_lookup_t type_index_lookup;
 | 
|---|
| 320 |   std::vector<unsigned int> elementnumbers;
 | 
|---|
| 321 |   elementnumbers.reserve(params.graph_elements.get().size());
 | 
|---|
| 322 |   size_t i = 0;
 | 
|---|
| 323 |   for (std::vector<const element *>::const_iterator iter = params.graph_elements.get().begin();
 | 
|---|
| 324 |       iter != params.graph_elements.get().end(); ++iter) {
 | 
|---|
| 325 |     elementnumbers.push_back((*iter)->getAtomicNumber());
 | 
|---|
| 326 |     LOG(3, "DEBUG: Adding [" << i << ", " << elementnumbers.back() << "] type to index lookup.");
 | 
|---|
| 327 |     type_index_lookup.left.insert(
 | 
|---|
| 328 |         std::make_pair(i++, Extractors::ParticleType_t(elementnumbers.back())));
 | 
|---|
| 329 |   }
 | 
|---|
| 330 |   if (i != boost::num_vertices(graph)) {
 | 
|---|
| 331 |     ELOG(1, "Not enough elements given.");
 | 
|---|
| 332 |     return Action::failure;
 | 
|---|
| 333 |   }
 | 
|---|
| 334 | 
 | 
|---|
| 335 |   /// 3a. get range of valence for each node
 | 
|---|
| 336 |   type_valence_map_t type_valence_map;
 | 
|---|
| 337 |   for (std::vector<const element *>::const_iterator iter = params.graph_elements.get().begin();
 | 
|---|
| 338 |       iter != params.graph_elements.get().end(); ++iter) {
 | 
|---|
| 339 |     if (type_valence_map.count((*iter)->getAtomicNumber()) == 0) {
 | 
|---|
| 340 |       LOG(3, "DEBUG: Adding [" << (*iter)->getAtomicNumber()
 | 
|---|
| 341 |           << ", " << (*iter)->getNoValenceOrbitals() << "] type to valence lookup.");
 | 
|---|
| 342 |       type_valence_map.insert( std::make_pair(
 | 
|---|
| 343 |           (*iter)->getAtomicNumber(),
 | 
|---|
| 344 |           (*iter)->getNoValenceOrbitals()) );
 | 
|---|
| 345 |     }
 | 
|---|
| 346 |   }
 | 
|---|
| 347 | 
 | 
|---|
| 348 |   // get map for retrieving vertex index for a given vertex
 | 
|---|
| 349 |   const BoostGraphCreator::name_map_t name_map = boost::get(boost::vertex_name, graph);
 | 
|---|
| 350 | 
 | 
|---|
| 351 |   /// 3b. generate range of degree for each edge
 | 
|---|
| 352 |   edge_valence_t edge_valence;
 | 
|---|
| 353 |   edges_t edges;
 | 
|---|
| 354 |   edge_index_t edge_index;
 | 
|---|
| 355 |   edges.reserve(boost::num_edges(graph));
 | 
|---|
| 356 |   {
 | 
|---|
| 357 |     size_t index = 0;
 | 
|---|
| 358 |     BoostGraphCreator::edge_iter i, end;
 | 
|---|
| 359 |     for (boost::tie(i, end) = boost::edges(graph); i != end; ++i) {
 | 
|---|
| 360 |       static int MAX_DEGREE = 3;
 | 
|---|
| 361 |       const BoostGraphCreator::Edge &e = *i;
 | 
|---|
| 362 |       /// 4a. place all edges into a vector
 | 
|---|
| 363 |       const BoostGraphCreator::Vertex u = source(e, graph);
 | 
|---|
| 364 |       const BoostGraphCreator::Vertex v = target(e, graph);
 | 
|---|
| 365 |       const edge_by_vertices_t edge_by_vertices = getEdgeByVertices(u,v);
 | 
|---|
| 366 |       edges.push_back(edge_by_vertices);
 | 
|---|
| 367 |       edge_index.insert( std::make_pair(
 | 
|---|
| 368 |           edge_by_vertices, index++) );
 | 
|---|
| 369 |       const int valence_u = getValenceForVertex(name_map, type_index_lookup, type_valence_map, u);
 | 
|---|
| 370 |       const int valence_v = getValenceForVertex(name_map, type_index_lookup, type_valence_map, v);
 | 
|---|
| 371 |       const int max_degree = std::max(1,
 | 
|---|
| 372 |           std::min(
 | 
|---|
| 373 |               std::min(valence_u, MAX_DEGREE),
 | 
|---|
| 374 |               std::min(valence_v, MAX_DEGREE)
 | 
|---|
| 375 |               ));
 | 
|---|
| 376 |       edge_valence.insert(
 | 
|---|
| 377 |           std::make_pair(edge_by_vertices, max_degree)
 | 
|---|
| 378 |           );
 | 
|---|
| 379 |       LOG(3, "DEBUG: Adding [" << std::make_pair(u,v) << ", "
 | 
|---|
| 380 |           << max_degree << "] edge to valence lookup.");
 | 
|---|
| 381 |     }
 | 
|---|
| 382 |   }
 | 
|---|
| 383 | 
 | 
|---|
| 384 |   /// 4. for each graph with type per node and degree per edge set, add hydrogen nodes and edges
 | 
|---|
| 385 |   // go through powerset over all edges
 | 
|---|
| 386 | 
 | 
|---|
| 387 |   /// 4b. have a recursive function that puts together the powerset over all max_degree
 | 
|---|
| 388 |   list_of_edge_degrees_t list_of_edge_degrees;
 | 
|---|
| 389 |   size_t no_combinations = 0;
 | 
|---|
| 390 |   {
 | 
|---|
| 391 |     size_t remaining_edges = edges.size();
 | 
|---|
| 392 |     size_t current_index = 0;
 | 
|---|
| 393 |     degrees_t current_degrees(edges.size(), 0);
 | 
|---|
| 394 |     no_combinations = powerset_edge_degree(
 | 
|---|
| 395 |         list_of_edge_degrees, current_degrees,
 | 
|---|
| 396 |         remaining_edges, current_index,
 | 
|---|
| 397 |         edges, edge_valence);
 | 
|---|
| 398 |   }
 | 
|---|
| 399 |   LOG(1, "INFO: Added " << no_combinations << " graph degree combinations.");
 | 
|---|
| 400 | 
 | 
|---|
| 401 |   /// 5. dissect graph into a list of subgraphs
 | 
|---|
| 402 |   typedef std::vector<BoostGraphCreator::UndirectedGraph> subgraphs_t;
 | 
|---|
| 403 |   subgraphs_t allsubgraphs;
 | 
|---|
| 404 |   const size_t nodes_in_graph = boost::num_vertices(graph);
 | 
|---|
| 405 |   LOG(2, "DEBUG: There are " << nodes_in_graph << " nodes in the graph.");
 | 
|---|
| 406 |   Extractors::set_of_nodes_t set_of_nodes;
 | 
|---|
| 407 |   for (size_t nodeid = 0; nodeid < nodes_in_graph; ++nodeid) {
 | 
|---|
| 408 |     const size_t &rootindex = nodeid;
 | 
|---|
| 409 |     LOG(2, "DEBUG: Current root index is " << rootindex);
 | 
|---|
| 410 | 
 | 
|---|
| 411 |     /// 5a. from node in graph (with this index) perform BFS till n-1 (#nodes in binding model)
 | 
|---|
| 412 |     std::vector<size_t> distances(boost::num_vertices(graph));
 | 
|---|
| 413 |     boost::breadth_first_search(
 | 
|---|
| 414 |         graph,
 | 
|---|
| 415 |         boost::vertex(rootindex, graph),
 | 
|---|
| 416 |         boost::visitor(Extractors::record_distance(&distances[0])));
 | 
|---|
| 417 |     LOG(3, "DEBUG: BFS discovered the following distances " << distances);
 | 
|---|
| 418 | 
 | 
|---|
| 419 |     /// 5b. and store all node in map with distance to root as key
 | 
|---|
| 420 |     Extractors::nodes_per_level_t nodes_per_level;
 | 
|---|
| 421 |     for (size_t i=0;i<distances.size();++i) {
 | 
|---|
| 422 |       nodes_per_level.insert( std::make_pair(Extractors::level_t(distances[i]), Extractors::node_t(i)) );
 | 
|---|
| 423 |     }
 | 
|---|
| 424 |     LOG(3, "DEBUG: Nodes per level are " << nodes_per_level);
 | 
|---|
| 425 | 
 | 
|---|
| 426 |     /// 5c. construct all possible induced connected subgraphs with this map (see fragmentation)
 | 
|---|
| 427 |     Extractors::nodes_t nodes;
 | 
|---|
| 428 |     // rootindex is always contained
 | 
|---|
| 429 |     nodes.insert( rootindex );
 | 
|---|
| 430 |     Extractors::level_t level = 0;
 | 
|---|
| 431 | 
 | 
|---|
| 432 |     for (size_t i=0;i<nodes_in_graph; ++i)
 | 
|---|
| 433 |       Extractors::generateAllInducedConnectedSubgraphs(
 | 
|---|
| 434 |           i, level, nodes, set_of_nodes, nodes_per_level, graph, distances, index_map);
 | 
|---|
| 435 |     LOG(2, "DEBUG: We have found " << set_of_nodes.size() << " unique induced, connected subgraphs.");
 | 
|---|
| 436 |   }
 | 
|---|
| 437 | 
 | 
|---|
| 438 |   // put all nodes in a list
 | 
|---|
| 439 |   Extractors::nodes_t nodes;
 | 
|---|
| 440 |   for (size_t nodeid = 0; nodeid < nodes_in_graph; ++nodeid) {
 | 
|---|
| 441 |     nodes.insert(nodeid);
 | 
|---|
| 442 |   }
 | 
|---|
| 443 | 
 | 
|---|
| 444 |   // 6. for every combination saturate fragments for lookup and energy contribution summation
 | 
|---|
| 445 |   size_t num_admissible = 0;
 | 
|---|
| 446 |   const HomologyContainer &container = World::getInstance().getHomologies();
 | 
|---|
| 447 |   for (list_of_edge_degrees_t::const_iterator degrees_iter = list_of_edge_degrees.begin();
 | 
|---|
| 448 |       degrees_iter != list_of_edge_degrees.end(); ++degrees_iter) {
 | 
|---|
| 449 |     const degrees_t ¤t_degrees = *degrees_iter;
 | 
|---|
| 450 |     LOG(2, "DEBUG: Current degree combination is " << current_degrees);
 | 
|---|
| 451 | 
 | 
|---|
| 452 |     // saturate subgraph: add to nodes, add to type_index_lookup, index_map
 | 
|---|
| 453 |       // saturate any leftover valence with hydrogen nodes
 | 
|---|
| 454 |     SaturatedGraph fullgraph =
 | 
|---|
| 455 |         saturateGraph(graph, name_map,
 | 
|---|
| 456 |             type_index_lookup, type_valence_map, edge_index,
 | 
|---|
| 457 |             edges, current_degrees);
 | 
|---|
| 458 | 
 | 
|---|
| 459 |     // don't add if one of the node's valence cannot be fulfilled
 | 
|---|
| 460 |     if (!fullgraph.ValencesAllFulfilled) {
 | 
|---|
| 461 |       LOG(2, "DEBUG: The degree combination " << current_degrees << " can NOT be fulfilled.");
 | 
|---|
| 462 |       continue;
 | 
|---|
| 463 |     } else {
 | 
|---|
| 464 |       ++num_admissible;
 | 
|---|
| 465 |       LOG(2, "DEBUG: The degree combination is admissable.");
 | 
|---|
| 466 |     }
 | 
|---|
| 467 | 
 | 
|---|
| 468 | 
 | 
|---|
| 469 |     // simply sum up the contributions from all fragments for the total energy
 | 
|---|
| 470 |     // of this particular combination, given the graph, the element for each node,
 | 
|---|
| 471 |     // and the bond degrees for each edge
 | 
|---|
| 472 |     double total_energy = 0.;
 | 
|---|
| 473 | 
 | 
|---|
| 474 |     // get fragment
 | 
|---|
| 475 |     for (Extractors::set_of_nodes_t::const_iterator nodesiter = set_of_nodes.begin();
 | 
|---|
| 476 |         nodesiter != set_of_nodes.end(); ++nodesiter) {
 | 
|---|
| 477 |       const Extractors::nodes_t ¤t_nodes = *nodesiter;
 | 
|---|
| 478 |       LOG(2, "DEBUG: Fragment nodes for graph are " << current_nodes);
 | 
|---|
| 479 | 
 | 
|---|
| 480 |       // create subgraph
 | 
|---|
| 481 |       edge_index_t subgraph_edge_index;
 | 
|---|
| 482 |       Extractors::UndirectedGraph newgraph = extractSubgraph(
 | 
|---|
| 483 |           graph, current_nodes, edge_index, subgraph_edge_index);
 | 
|---|
| 484 | 
 | 
|---|
| 485 |       const BoostGraphCreator::name_map_t new_name_map =
 | 
|---|
| 486 |           boost::get(boost::vertex_name, newgraph);
 | 
|---|
| 487 |       {
 | 
|---|
| 488 |         BoostGraphCreator::vertex_iter viter, vend;
 | 
|---|
| 489 |         boost::tie(viter, vend) = boost::vertices(newgraph);
 | 
|---|
| 490 |         for (; viter != vend; ++viter) {
 | 
|---|
| 491 |           const atomId_t vindex = boost::get(new_name_map, *viter);
 | 
|---|
| 492 |           LOG(4, "DEBUG: Copied node in fragment has index " << vindex);
 | 
|---|
| 493 |         }
 | 
|---|
| 494 |       }
 | 
|---|
| 495 | 
 | 
|---|
| 496 |       // saturate subgraph: add to nodes, add to type_index_lookup, index_map
 | 
|---|
| 497 |         // saturate any leftover valence with hydrogen nodes
 | 
|---|
| 498 |       SaturatedGraph fragmentgraph =
 | 
|---|
| 499 |           saturateGraph(newgraph, new_name_map,
 | 
|---|
| 500 |               type_index_lookup, type_valence_map, subgraph_edge_index,
 | 
|---|
| 501 |               edges, current_degrees);
 | 
|---|
| 502 | 
 | 
|---|
| 503 |       // don't add if one of the node's valence cannot be fulfilled
 | 
|---|
| 504 |       ASSERT(fragmentgraph.ValencesAllFulfilled,
 | 
|---|
| 505 |           "GraphChemicalSpaceEvaluatorAction::performCall() - encountered subgraph whose valences cannot be fulfilled.");
 | 
|---|
| 506 | 
 | 
|---|
| 507 |       const Extractors::index_map_t subgraph_index_map =
 | 
|---|
| 508 |           boost::get(boost::vertex_index, fragmentgraph.saturated_graph);
 | 
|---|
| 509 |       // const BoostGraphCreator::name_map_t subgraph_name_map =
 | 
|---|
| 510 |       //    boost::get(boost::vertex_name, fragmentgraph.saturated_graph);
 | 
|---|
| 511 | 
 | 
|---|
| 512 |       Extractors::nodes_t saturated_nodes;
 | 
|---|
| 513 |       BoostGraphCreator::vertex_iter viter, vend;
 | 
|---|
| 514 |       boost::tie(viter, vend) = boost::vertices(fragmentgraph.saturated_graph);
 | 
|---|
| 515 |       for (; viter != vend; ++viter) {
 | 
|---|
| 516 |         const Extractors::node_t vindex = boost::get(subgraph_index_map, *viter);
 | 
|---|
| 517 |         LOG(4, "DEBUG: Adding fragment node index " << vindex << " of type "
 | 
|---|
| 518 |             << Extractors::getParticleTypeToNode(fragmentgraph.type_index_lookup, vindex)
 | 
|---|
| 519 |             << " to saturated fragment nodes list.");
 | 
|---|
| 520 |         saturated_nodes.insert(vindex);
 | 
|---|
| 521 |       }
 | 
|---|
| 522 |       /// convert its nodes into a HomologyGraph using the graph and elements (types for index)
 | 
|---|
| 523 |       const HomologyGraph nodes_graph =
 | 
|---|
| 524 |           Extractors::createHomologyGraphFromNodes(
 | 
|---|
| 525 |               saturated_nodes,
 | 
|---|
| 526 |               fragmentgraph.type_index_lookup,
 | 
|---|
| 527 |               fragmentgraph.saturated_graph,
 | 
|---|
| 528 |               subgraph_index_map);
 | 
|---|
| 529 |       LOG(2, "DEBUG: Looking for fragment graphs " << nodes_graph << " in homology container.");
 | 
|---|
| 530 | 
 | 
|---|
| 531 |       /// Query HomologyContainer for this HomologyGraph.
 | 
|---|
| 532 |       HomologyContainer::range_t range = container.getHomologousGraphs(nodes_graph);
 | 
|---|
| 533 | 
 | 
|---|
| 534 |       if (range.first == range.second) {
 | 
|---|
| 535 |         // range is empty, the fragment is unknown
 | 
|---|
| 536 |         ELOG(1, "Cannot find fragment graph " << nodes_graph << " to graph " << elementnumbers);
 | 
|---|
| 537 |       } else {
 | 
|---|
| 538 |         // list lowest energy
 | 
|---|
| 539 |         const HomologyContainer::const_iterator lowest_contribution_graph =
 | 
|---|
| 540 |             std::min_element(range.first, range.second, compareEnergyContribution);
 | 
|---|
| 541 |         const HomologyContainer::const_iterator highest_contribution_graph =
 | 
|---|
| 542 |             std::max_element(range.first, range.second, compareEnergyContribution);
 | 
|---|
| 543 |         LOG(2, "INFO: Fragment graph " << nodes_graph << " has energy contributions from "
 | 
|---|
| 544 |             << lowest_contribution_graph->second.energy << " Ht till "
 | 
|---|
| 545 |             << highest_contribution_graph->second.energy << " Ht, picking lowest.");
 | 
|---|
| 546 |         total_energy += lowest_contribution_graph->second.energy;
 | 
|---|
| 547 |       }
 | 
|---|
| 548 |     }
 | 
|---|
| 549 |     LOG(1, "The graph with degrees " << current_degrees << " has a total BOSSANOVA energy of " << total_energy);
 | 
|---|
| 550 |   }
 | 
|---|
| 551 |   LOG(1, "There have been " << num_admissible << " admissible degree combinations for the given graph.");
 | 
|---|
| 552 | 
 | 
|---|
| 553 |   return Action::success;
 | 
|---|
| 554 | }
 | 
|---|
| 555 | 
 | 
|---|
| 556 | ActionState::ptr GraphChemicalSpaceEvaluatorAction::performUndo(ActionState::ptr _state) {
 | 
|---|
| 557 |   return Action::success;
 | 
|---|
| 558 | }
 | 
|---|
| 559 | 
 | 
|---|
| 560 | ActionState::ptr GraphChemicalSpaceEvaluatorAction::performRedo(ActionState::ptr _state){
 | 
|---|
| 561 |   return Action::success;
 | 
|---|
| 562 | }
 | 
|---|
| 563 | 
 | 
|---|
| 564 | bool GraphChemicalSpaceEvaluatorAction::canUndo() {
 | 
|---|
| 565 |   return true;
 | 
|---|
| 566 | }
 | 
|---|
| 567 | 
 | 
|---|
| 568 | bool GraphChemicalSpaceEvaluatorAction::shouldUndo() {
 | 
|---|
| 569 |   return true;
 | 
|---|
| 570 | }
 | 
|---|
| 571 | /** =========== end of function ====================== */
 | 
|---|