| [f5ea10] | 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>
 | 
|---|
| [99c705] | 38 | #include <list>
 | 
|---|
 | 39 | #include <map>
 | 
|---|
| [f5ea10] | 40 | #include <string>
 | 
|---|
 | 41 | #include <vector>
 | 
|---|
 | 42 | 
 | 
|---|
| [99c705] | 43 | #include <boost/graph/copy.hpp>
 | 
|---|
 | 44 | 
 | 
|---|
| [f5ea10] | 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 ====================== */
 | 
|---|
| [99c705] | 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);
 | 
|---|
| [1be100] | 144 |   const BoostGraphCreator::Vertex v = target(e, graph);
 | 
|---|
| [99c705] | 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 | 
 | 
|---|
| [d951a5] | 193 |   LOG(2, "DEBUG: Saturating graph with " << boost::num_vertices(newgraph.saturated_graph) << " nodes.");
 | 
|---|
| [99c705] | 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;
 | 
|---|
| [1be100] | 206 |       const edge_by_vertices_t edge_by_vertices = getEdgeVerticesByEdge(e, graph);
 | 
|---|
| [99c705] | 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);
 | 
|---|
| [1be100] | 231 |       ASSERT( edge_inserter.second,
 | 
|---|
 | 232 |           "Failed to insert hydrogen into saturation graph.");
 | 
|---|
| [99c705] | 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,
 | 
|---|
| [1be100] | 244 |     const Extractors::nodes_t &nodes,
 | 
|---|
 | 245 |     const edge_index_t &_edge_index,
 | 
|---|
 | 246 |     edge_index_t &_subgraph_edge_index)
 | 
|---|
| [99c705] | 247 | {
 | 
|---|
 | 248 |   BoostGraphCreator::UndirectedGraph subgraph;
 | 
|---|
 | 249 |   const Extractors::index_map_t index_map = boost::get(boost::vertex_index, _graph);
 | 
|---|
| [1be100] | 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;
 | 
|---|
| [99c705] | 252 | 
 | 
|---|
 | 253 |   // add nodes
 | 
|---|
 | 254 |   BoostGraphCreator::vertex_iter viter, vend;
 | 
|---|
 | 255 |   for (boost::tie(viter, vend) = boost::vertices(_graph); viter != vend; ++viter) {
 | 
|---|
| [1be100] | 256 |     const Extractors::node_t nodeid = boost::get(index_map, *viter);
 | 
|---|
| [99c705] | 257 |     if (nodes.find(nodeid) != nodes.end()) {
 | 
|---|
| [1be100] | 258 |       const BoostGraphCreator::Vertex vnew = boost::add_vertex(*viter, subgraph);
 | 
|---|
 | 259 |       graph_index_to_subgraph_vertex.insert( std::make_pair(nodeid, vnew) );
 | 
|---|
| [99c705] | 260 |       LOG(4, "DEBUG: Adding node " << *viter << " to subgraph.");
 | 
|---|
 | 261 |     }
 | 
|---|
 | 262 |   }
 | 
|---|
| [1be100] | 263 |   const Extractors::index_map_t subgraph_index_map = boost::get(boost::vertex_index, subgraph);
 | 
|---|
| [99c705] | 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));
 | 
|---|
| [1be100] | 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]);
 | 
|---|
| [99c705] | 277 |       if ((nodes.find(sourceindex) != nodes.end()) && (nodes.find(targetindex) != nodes.end())
 | 
|---|
 | 278 |           && (sourceindex < targetindex)) {
 | 
|---|
| [1be100] | 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 << ".");
 | 
|---|
| [99c705] | 294 |       }
 | 
|---|
 | 295 |     }
 | 
|---|
 | 296 |   }
 | 
|---|
 | 297 | 
 | 
|---|
 | 298 |   return subgraph;
 | 
|---|
 | 299 | }
 | 
|---|
 | 300 | 
 | 
|---|
| [f5ea10] | 301 | ActionState::ptr GraphChemicalSpaceEvaluatorAction::performCall() {
 | 
|---|
| [99c705] | 302 |   /// 1. create boost::graph from graph6 string
 | 
|---|
| [f5ea10] | 303 |   BoostGraphCreator BGC;
 | 
|---|
 | 304 |   BGC.createFromGraph6String(params.graph_string.get());
 | 
|---|
 | 305 | 
 | 
|---|
 | 306 |   BoostGraphCreator::UndirectedGraph graph = BGC.get();
 | 
|---|
 | 307 | 
 | 
|---|
 | 308 |   Extractors::index_map_t index_map = boost::get(boost::vertex_index, graph);
 | 
|---|
| [99c705] | 309 |   LOG(1, "INFO: We have " << boost::num_vertices(graph) << " nodes in the fragment graph.");
 | 
|---|
| [f5ea10] | 310 | 
 | 
|---|
| [99c705] | 311 |   /// 2. associate type with a node index
 | 
|---|
| [f5ea10] | 312 |   Extractors::type_index_lookup_t type_index_lookup;
 | 
|---|
 | 313 |   std::vector<unsigned int> elementnumbers;
 | 
|---|
 | 314 |   elementnumbers.reserve(params.graph_elements.get().size());
 | 
|---|
| [99c705] | 315 |   size_t i = 0;
 | 
|---|
| [f5ea10] | 316 |   for (std::vector<const element *>::const_iterator iter = params.graph_elements.get().begin();
 | 
|---|
 | 317 |       iter != params.graph_elements.get().end(); ++iter) {
 | 
|---|
 | 318 |     elementnumbers.push_back((*iter)->getAtomicNumber());
 | 
|---|
| [99c705] | 319 |     LOG(3, "DEBUG: Adding [" << i << ", " << elementnumbers.back() << "] type to index lookup.");
 | 
|---|
| [f5ea10] | 320 |     type_index_lookup.left.insert(
 | 
|---|
| [99c705] | 321 |         std::make_pair(i++, Extractors::ParticleType_t(elementnumbers.back())));
 | 
|---|
 | 322 |   }
 | 
|---|
 | 323 |   if (i != boost::num_vertices(graph)) {
 | 
|---|
 | 324 |     ELOG(1, "Not enough elements given.");
 | 
|---|
 | 325 |     return Action::failure;
 | 
|---|
| [f5ea10] | 326 |   }
 | 
|---|
 | 327 | 
 | 
|---|
| [99c705] | 328 |   /// 3a. get range of valence for each node
 | 
|---|
 | 329 |   type_valence_map_t type_valence_map;
 | 
|---|
 | 330 |   for (std::vector<const element *>::const_iterator iter = params.graph_elements.get().begin();
 | 
|---|
 | 331 |       iter != params.graph_elements.get().end(); ++iter) {
 | 
|---|
 | 332 |     if (type_valence_map.count((*iter)->getAtomicNumber()) == 0) {
 | 
|---|
 | 333 |       LOG(3, "DEBUG: Adding [" << (*iter)->getAtomicNumber()
 | 
|---|
 | 334 |           << ", " << (*iter)->getNoValenceOrbitals() << "] type to valence lookup.");
 | 
|---|
 | 335 |       type_valence_map.insert( std::make_pair(
 | 
|---|
 | 336 |           (*iter)->getAtomicNumber(),
 | 
|---|
 | 337 |           (*iter)->getNoValenceOrbitals()) );
 | 
|---|
 | 338 |     }
 | 
|---|
 | 339 |   }
 | 
|---|
 | 340 | 
 | 
|---|
 | 341 |   // get map for retrieving vertex index for a given vertex
 | 
|---|
 | 342 |   const BoostGraphCreator::name_map_t name_map = boost::get(boost::vertex_name, graph);
 | 
|---|
 | 343 | 
 | 
|---|
 | 344 |   /// 3b. generate range of degree for each edge
 | 
|---|
 | 345 |   edge_valence_t edge_valence;
 | 
|---|
 | 346 |   edges_t edges;
 | 
|---|
 | 347 |   edge_index_t edge_index;
 | 
|---|
 | 348 |   edges.reserve(boost::num_edges(graph));
 | 
|---|
 | 349 |   {
 | 
|---|
 | 350 |     size_t index = 0;
 | 
|---|
 | 351 |     BoostGraphCreator::edge_iter i, end;
 | 
|---|
 | 352 |     for (boost::tie(i, end) = boost::edges(graph); i != end; ++i) {
 | 
|---|
 | 353 |       static int MAX_DEGREE = 3;
 | 
|---|
 | 354 |       const BoostGraphCreator::Edge &e = *i;
 | 
|---|
 | 355 |       /// 4a. place all edges into a vector
 | 
|---|
 | 356 |       const BoostGraphCreator::Vertex u = source(e, graph);
 | 
|---|
 | 357 |       const BoostGraphCreator::Vertex v = target(e, graph);
 | 
|---|
 | 358 |       const edge_by_vertices_t edge_by_vertices = getEdgeByVertices(u,v);
 | 
|---|
 | 359 |       edges.push_back(edge_by_vertices);
 | 
|---|
 | 360 |       edge_index.insert( std::make_pair(
 | 
|---|
 | 361 |           edge_by_vertices, index++) );
 | 
|---|
 | 362 |       const int valence_u = getValenceForVertex(name_map, type_index_lookup, type_valence_map, u);
 | 
|---|
 | 363 |       const int valence_v = getValenceForVertex(name_map, type_index_lookup, type_valence_map, v);
 | 
|---|
 | 364 |       const int max_degree = std::max(1,
 | 
|---|
 | 365 |           std::min(
 | 
|---|
 | 366 |               std::min(valence_u, MAX_DEGREE),
 | 
|---|
 | 367 |               std::min(valence_v, MAX_DEGREE)
 | 
|---|
 | 368 |               ));
 | 
|---|
 | 369 |       edge_valence.insert(
 | 
|---|
 | 370 |           std::make_pair(edge_by_vertices, max_degree)
 | 
|---|
 | 371 |           );
 | 
|---|
 | 372 |       LOG(3, "DEBUG: Adding [" << std::make_pair(u,v) << ", "
 | 
|---|
 | 373 |           << max_degree << "] edge to valence lookup.");
 | 
|---|
 | 374 |     }
 | 
|---|
 | 375 |   }
 | 
|---|
 | 376 | 
 | 
|---|
 | 377 |   /// 4. for each graph with type per node and degree per edge set, add hydrogen nodes and edges
 | 
|---|
 | 378 |   // go through powerset over all edges
 | 
|---|
 | 379 | 
 | 
|---|
 | 380 |   /// 4b. have a recursive function that puts together the powerset over all max_degree
 | 
|---|
 | 381 |   list_of_edge_degrees_t list_of_edge_degrees;
 | 
|---|
 | 382 |   size_t no_combinations = 0;
 | 
|---|
 | 383 |   {
 | 
|---|
 | 384 |     size_t remaining_edges = edges.size();
 | 
|---|
 | 385 |     size_t current_index = 0;
 | 
|---|
 | 386 |     degrees_t current_degrees(edges.size(), 0);
 | 
|---|
 | 387 |     no_combinations = powerset_edge_degree(
 | 
|---|
 | 388 |         list_of_edge_degrees, current_degrees,
 | 
|---|
 | 389 |         remaining_edges, current_index,
 | 
|---|
 | 390 |         edges, edge_valence);
 | 
|---|
 | 391 |   }
 | 
|---|
 | 392 |   LOG(1, "INFO: Added " << no_combinations << " graph degree combinations.");
 | 
|---|
 | 393 | 
 | 
|---|
 | 394 |   /// 5. dissect graph into a list of subgraphs
 | 
|---|
 | 395 |   typedef std::vector<BoostGraphCreator::UndirectedGraph> subgraphs_t;
 | 
|---|
 | 396 |   subgraphs_t allsubgraphs;
 | 
|---|
| [f5ea10] | 397 |   const size_t nodes_in_graph = boost::num_vertices(graph);
 | 
|---|
 | 398 |   LOG(2, "DEBUG: There are " << nodes_in_graph << " nodes in the graph.");
 | 
|---|
 | 399 |   Extractors::set_of_nodes_t set_of_nodes;
 | 
|---|
 | 400 |   for (size_t nodeid = 0; nodeid < nodes_in_graph; ++nodeid) {
 | 
|---|
 | 401 |     const size_t &rootindex = nodeid;
 | 
|---|
 | 402 |     LOG(2, "DEBUG: Current root index is " << rootindex);
 | 
|---|
 | 403 | 
 | 
|---|
 | 404 |     /// 5a. from node in graph (with this index) perform BFS till n-1 (#nodes in binding model)
 | 
|---|
 | 405 |     std::vector<size_t> distances(boost::num_vertices(graph));
 | 
|---|
 | 406 |     boost::breadth_first_search(
 | 
|---|
 | 407 |         graph,
 | 
|---|
 | 408 |         boost::vertex(rootindex, graph),
 | 
|---|
 | 409 |         boost::visitor(Extractors::record_distance(&distances[0])));
 | 
|---|
 | 410 |     LOG(3, "DEBUG: BFS discovered the following distances " << distances);
 | 
|---|
 | 411 | 
 | 
|---|
 | 412 |     /// 5b. and store all node in map with distance to root as key
 | 
|---|
 | 413 |     Extractors::nodes_per_level_t nodes_per_level;
 | 
|---|
 | 414 |     for (size_t i=0;i<distances.size();++i) {
 | 
|---|
 | 415 |       nodes_per_level.insert( std::make_pair(Extractors::level_t(distances[i]), Extractors::node_t(i)) );
 | 
|---|
 | 416 |     }
 | 
|---|
 | 417 |     LOG(3, "DEBUG: Nodes per level are " << nodes_per_level);
 | 
|---|
 | 418 | 
 | 
|---|
| [99c705] | 419 |     /// 5c. construct all possible induced connected subgraphs with this map (see fragmentation)
 | 
|---|
| [f5ea10] | 420 |     Extractors::nodes_t nodes;
 | 
|---|
 | 421 |     // rootindex is always contained
 | 
|---|
 | 422 |     nodes.insert( rootindex );
 | 
|---|
 | 423 |     Extractors::level_t level = 0;
 | 
|---|
 | 424 | 
 | 
|---|
| [99c705] | 425 |     for (size_t i=0;i<nodes_in_graph; ++i)
 | 
|---|
 | 426 |       Extractors::generateAllInducedConnectedSubgraphs(
 | 
|---|
 | 427 |           i, level, nodes, set_of_nodes, nodes_per_level, graph, distances, index_map);
 | 
|---|
| [f5ea10] | 428 |     LOG(2, "DEBUG: We have found " << set_of_nodes.size() << " unique induced, connected subgraphs.");
 | 
|---|
 | 429 |   }
 | 
|---|
 | 430 | 
 | 
|---|
| [d951a5] | 431 |   // put all nodes in a list
 | 
|---|
 | 432 |   Extractors::nodes_t nodes;
 | 
|---|
 | 433 |   for (size_t nodeid = 0; nodeid < nodes_in_graph; ++nodeid) {
 | 
|---|
 | 434 |     nodes.insert(nodeid);
 | 
|---|
 | 435 |   }
 | 
|---|
 | 436 | 
 | 
|---|
| [99c705] | 437 |   // 6. for every combination saturate fragments for lookup and energy contribution summation
 | 
|---|
| [1be100] | 438 |   size_t num_admissible = 0;
 | 
|---|
| [f5ea10] | 439 |   const HomologyContainer &container = World::getInstance().getHomologies();
 | 
|---|
| [99c705] | 440 |   for (list_of_edge_degrees_t::const_iterator degrees_iter = list_of_edge_degrees.begin();
 | 
|---|
 | 441 |       degrees_iter != list_of_edge_degrees.end(); ++degrees_iter) {
 | 
|---|
 | 442 |     const degrees_t ¤t_degrees = *degrees_iter;
 | 
|---|
| [d951a5] | 443 |     LOG(2, "DEBUG: Current degree combination is " << current_degrees);
 | 
|---|
 | 444 | 
 | 
|---|
 | 445 |     // saturate subgraph: add to nodes, add to type_index_lookup, index_map
 | 
|---|
 | 446 |       // saturate any leftover valence with hydrogen nodes
 | 
|---|
 | 447 |     SaturatedGraph fullgraph =
 | 
|---|
 | 448 |         saturateGraph(graph, name_map,
 | 
|---|
 | 449 |             type_index_lookup, type_valence_map, edge_index,
 | 
|---|
 | 450 |             edges, current_degrees);
 | 
|---|
 | 451 | 
 | 
|---|
 | 452 |     // don't add if one of the node's valence cannot be fulfilled
 | 
|---|
 | 453 |     if (!fullgraph.ValencesAllFulfilled) {
 | 
|---|
 | 454 |       LOG(2, "DEBUG: The degree combination " << current_degrees << " can NOT be fulfilled.");
 | 
|---|
 | 455 |       continue;
 | 
|---|
 | 456 |     } else {
 | 
|---|
| [1be100] | 457 |       ++num_admissible;
 | 
|---|
| [d951a5] | 458 |       LOG(2, "DEBUG: The degree combination is admissable.");
 | 
|---|
 | 459 |     }
 | 
|---|
 | 460 | 
 | 
|---|
| [99c705] | 461 | 
 | 
|---|
 | 462 |     // simply sum up the contributions from all fragments for the total energy
 | 
|---|
 | 463 |     // of this particular combination, given the graph, the element for each node,
 | 
|---|
 | 464 |     // and the bond degrees for each edge
 | 
|---|
 | 465 |     double total_energy = 0.;
 | 
|---|
 | 466 | 
 | 
|---|
 | 467 |     // get fragment
 | 
|---|
 | 468 |     for (Extractors::set_of_nodes_t::const_iterator nodesiter = set_of_nodes.begin();
 | 
|---|
 | 469 |         nodesiter != set_of_nodes.end(); ++nodesiter) {
 | 
|---|
 | 470 |       const Extractors::nodes_t ¤t_nodes = *nodesiter;
 | 
|---|
| [d951a5] | 471 |       LOG(2, "DEBUG: Fragment nodes for graph are " << current_nodes);
 | 
|---|
| [99c705] | 472 | 
 | 
|---|
 | 473 |       // create subgraph
 | 
|---|
| [1be100] | 474 |       edge_index_t subgraph_edge_index;
 | 
|---|
 | 475 |       Extractors::UndirectedGraph newgraph = extractSubgraph(
 | 
|---|
 | 476 |           graph, current_nodes, edge_index, subgraph_edge_index);
 | 
|---|
| [99c705] | 477 | 
 | 
|---|
 | 478 |       const BoostGraphCreator::name_map_t new_name_map =
 | 
|---|
 | 479 |           boost::get(boost::vertex_name, newgraph);
 | 
|---|
 | 480 |       {
 | 
|---|
 | 481 |         BoostGraphCreator::vertex_iter viter, vend;
 | 
|---|
 | 482 |         boost::tie(viter, vend) = boost::vertices(newgraph);
 | 
|---|
 | 483 |         for (; viter != vend; ++viter) {
 | 
|---|
 | 484 |           const atomId_t vindex = boost::get(new_name_map, *viter);
 | 
|---|
| [d951a5] | 485 |           LOG(4, "DEBUG: Copied node in fragment has index " << vindex);
 | 
|---|
| [99c705] | 486 |         }
 | 
|---|
 | 487 |       }
 | 
|---|
 | 488 | 
 | 
|---|
 | 489 |       // saturate subgraph: add to nodes, add to type_index_lookup, index_map
 | 
|---|
 | 490 |         // saturate any leftover valence with hydrogen nodes
 | 
|---|
| [d951a5] | 491 |       SaturatedGraph fragmentgraph =
 | 
|---|
 | 492 |           saturateGraph(newgraph, new_name_map,
 | 
|---|
| [1be100] | 493 |               type_index_lookup, type_valence_map, subgraph_edge_index,
 | 
|---|
| [99c705] | 494 |               edges, current_degrees);
 | 
|---|
 | 495 | 
 | 
|---|
 | 496 |       // don't add if one of the node's valence cannot be fulfilled
 | 
|---|
| [d951a5] | 497 |       ASSERT(fragmentgraph.ValencesAllFulfilled,
 | 
|---|
 | 498 |           "GraphChemicalSpaceEvaluatorAction::performCall() - encountered subgraph whose valences cannot be fulfilled.");
 | 
|---|
| [99c705] | 499 | 
 | 
|---|
 | 500 |       const Extractors::index_map_t subgraph_index_map =
 | 
|---|
| [d951a5] | 501 |           boost::get(boost::vertex_index, fragmentgraph.saturated_graph);
 | 
|---|
 | 502 |       // const BoostGraphCreator::name_map_t subgraph_name_map =
 | 
|---|
 | 503 |       //    boost::get(boost::vertex_name, fragmentgraph.saturated_graph);
 | 
|---|
| [99c705] | 504 | 
 | 
|---|
 | 505 |       Extractors::nodes_t saturated_nodes;
 | 
|---|
 | 506 |       BoostGraphCreator::vertex_iter viter, vend;
 | 
|---|
| [d951a5] | 507 |       boost::tie(viter, vend) = boost::vertices(fragmentgraph.saturated_graph);
 | 
|---|
| [99c705] | 508 |       for (; viter != vend; ++viter) {
 | 
|---|
 | 509 |         const Extractors::node_t vindex = boost::get(subgraph_index_map, *viter);
 | 
|---|
| [d951a5] | 510 |         LOG(4, "DEBUG: Adding fragment node index " << vindex << " of type "
 | 
|---|
 | 511 |             << Extractors::getParticleTypeToNode(fragmentgraph.type_index_lookup, vindex)
 | 
|---|
 | 512 |             << " to saturated fragment nodes list.");
 | 
|---|
| [99c705] | 513 |         saturated_nodes.insert(vindex);
 | 
|---|
 | 514 |       }
 | 
|---|
 | 515 |       /// convert its nodes into a HomologyGraph using the graph and elements (types for index)
 | 
|---|
 | 516 |       const HomologyGraph nodes_graph =
 | 
|---|
 | 517 |           Extractors::createHomologyGraphFromNodes(
 | 
|---|
 | 518 |               saturated_nodes,
 | 
|---|
| [d951a5] | 519 |               fragmentgraph.type_index_lookup,
 | 
|---|
 | 520 |               fragmentgraph.saturated_graph,
 | 
|---|
| [99c705] | 521 |               subgraph_index_map);
 | 
|---|
| [d951a5] | 522 |       LOG(2, "DEBUG: Looking for fragment graphs " << nodes_graph << " in homology container.");
 | 
|---|
| [99c705] | 523 | 
 | 
|---|
 | 524 |       /// Query HomologyContainer for this HomologyGraph.
 | 
|---|
 | 525 |       HomologyContainer::range_t range = container.getHomologousGraphs(nodes_graph);
 | 
|---|
 | 526 | 
 | 
|---|
 | 527 |       if (range.first == range.second) {
 | 
|---|
 | 528 |         // range is empty, the fragment is unknown
 | 
|---|
 | 529 |         ELOG(1, "Cannot find fragment graph " << nodes_graph << " to graph " << elementnumbers);
 | 
|---|
 | 530 |       } else {
 | 
|---|
| [9171d8] | 531 |         // list lowest energy
 | 
|---|
 | 532 |         const HomologyContainer::const_iterator lowest_contribution_graph =
 | 
|---|
| [ff2c52] | 533 |             std::min_element(range.first, range.second, HomologyContainer::compareEnergyContribution);
 | 
|---|
| [9171d8] | 534 |         const HomologyContainer::const_iterator highest_contribution_graph =
 | 
|---|
| [ff2c52] | 535 |             std::max_element(range.first, range.second, HomologyContainer::compareEnergyContribution);
 | 
|---|
| [9171d8] | 536 |         LOG(2, "INFO: Fragment graph " << nodes_graph << " has energy contributions from "
 | 
|---|
| [564f17] | 537 |             << lowest_contribution_graph->second.contribution << " Ht till "
 | 
|---|
 | 538 |             << highest_contribution_graph->second.contribution << " Ht, picking lowest.");
 | 
|---|
 | 539 |         total_energy += lowest_contribution_graph->second.contribution;
 | 
|---|
| [99c705] | 540 |       }
 | 
|---|
| [f5ea10] | 541 |     }
 | 
|---|
| [d951a5] | 542 |     LOG(1, "The graph with degrees " << current_degrees << " has a total BOSSANOVA energy of " << total_energy);
 | 
|---|
| [f5ea10] | 543 |   }
 | 
|---|
| [1be100] | 544 |   LOG(1, "There have been " << num_admissible << " admissible degree combinations for the given graph.");
 | 
|---|
| [f5ea10] | 545 | 
 | 
|---|
 | 546 |   return Action::success;
 | 
|---|
 | 547 | }
 | 
|---|
 | 548 | 
 | 
|---|
 | 549 | ActionState::ptr GraphChemicalSpaceEvaluatorAction::performUndo(ActionState::ptr _state) {
 | 
|---|
 | 550 |   return Action::success;
 | 
|---|
 | 551 | }
 | 
|---|
 | 552 | 
 | 
|---|
 | 553 | ActionState::ptr GraphChemicalSpaceEvaluatorAction::performRedo(ActionState::ptr _state){
 | 
|---|
 | 554 |   return Action::success;
 | 
|---|
 | 555 | }
 | 
|---|
 | 556 | 
 | 
|---|
 | 557 | bool GraphChemicalSpaceEvaluatorAction::canUndo() {
 | 
|---|
 | 558 |   return true;
 | 
|---|
 | 559 | }
 | 
|---|
 | 560 | 
 | 
|---|
 | 561 | bool GraphChemicalSpaceEvaluatorAction::shouldUndo() {
 | 
|---|
 | 562 |   return true;
 | 
|---|
 | 563 | }
 | 
|---|
 | 564 | /** =========== end of function ====================== */
 | 
|---|