/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2017 Frederik Heber. All rights reserved. * * * This file is part of MoleCuilder. * * MoleCuilder is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 2 of the License, or * (at your option) any later version. * * MoleCuilder is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with MoleCuilder. If not, see . */ /* * ChemicalSpaceEvaluator.cpp * * Created on: Sep 26, 2017 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif //#include "CodePatterns/MemDebug.hpp" #include #include #include #include #include #include #include "CodePatterns/Log.hpp" #include "Actions/GraphAction/ChemicalSpaceEvaluatorAction.hpp" #include "Element/element.hpp" #include "Graph/BoostGraphCreator.hpp" #include "Fragmentation/Homology/HomologyContainer.hpp" #include "FunctionApproximation/Extractors.hpp" #include "World.hpp" #include using namespace MoleCuilder; // and construct the stuff #include "ChemicalSpaceEvaluatorAction.def" #include "Action_impl_pre.hpp" /** =========== define the function ====================== */ template 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 edge_by_vertices_t; typedef std::map edge_valence_t; typedef std::vector edges_t; typedef std::vector 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 ¤t_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, 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 = target(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 static SaturatedGraph saturateGraph( const graph_type &graph, 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 °rees ) { 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); } } } LOG(2, "DEBUG: Saturating graph with " << boost::num_vertices(newgraph.saturated_graph) << " nodes."); 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::out_edge_iterator i, end; for (boost::tie(i, end) = boost::out_edges(v, graph); i != end; ++i) { const BoostGraphCreator::Edge &e = *i; const edge_by_vertices_t edge_by_vertices = getEdgeVerticesByEdge(e, graph); LOG(3, "DEBUG: Current edge is (" << edge_by_vertices << ")"); const size_t index = getElementFromMap( 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 edge_inserter = boost::add_edge(v, h, newgraph.saturated_graph); ASSERT( edge_inserter.second, "Failed to insert hydrogen into saturation graph."); } LOG(2, "DEBUG: Added " << (max_valence - total_valence) << " hydrogens to vertex #" << vindex); } return newgraph; } template BoostGraphCreator::UndirectedGraph extractSubgraph( const graph_t &_graph, const Extractors::nodes_t &nodes, const edge_index_t &_edge_index, edge_index_t &_subgraph_edge_index) { BoostGraphCreator::UndirectedGraph subgraph; const Extractors::index_map_t index_map = boost::get(boost::vertex_index, _graph); typedef std::map graph_index_to_subgraph_vertex_t; graph_index_to_subgraph_vertex_t graph_index_to_subgraph_vertex; // add nodes BoostGraphCreator::vertex_iter viter, vend; for (boost::tie(viter, vend) = boost::vertices(_graph); viter != vend; ++viter) { const Extractors::node_t nodeid = boost::get(index_map, *viter); if (nodes.find(nodeid) != nodes.end()) { const BoostGraphCreator::Vertex vnew = boost::add_vertex(*viter, subgraph); graph_index_to_subgraph_vertex.insert( std::make_pair(nodeid, vnew) ); LOG(4, "DEBUG: Adding node " << *viter << " to subgraph."); } } const Extractors::index_map_t subgraph_index_map = boost::get(boost::vertex_index, 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)); // we need to translate the vertex index from graph to subgraph const Extractors::node_t subsourceindex = boost::get( subgraph_index_map, graph_index_to_subgraph_vertex[sourceindex]); const Extractors::node_t subtargetindex = boost::get( subgraph_index_map, graph_index_to_subgraph_vertex[targetindex]); if ((nodes.find(sourceindex) != nodes.end()) && (nodes.find(targetindex) != nodes.end()) && (sourceindex < targetindex)) { // and we need to translate the edge index from the graph to the subgraph const std::pair newedgeinserter = boost::add_edge(subsourceindex, subtargetindex, subgraph); ASSERT( newedgeinserter.second, "extractSubgraph() - could not insert edge "+toString(subsourceindex)+"<->" +toString(subtargetindex)+"."); const edge_by_vertices_t edge_by_vertices = getEdgeVerticesByEdge(*i, _graph); const edge_index_t::const_iterator edgeiter = _edge_index.find(edge_by_vertices); ASSERT( edgeiter != _edge_index.end(), "extractSubgraph() - could not find edge "+toString(edge_by_vertices)+" in edge_index map." ); const edge_by_vertices_t subgraph_edge_by_vertices = getEdgeVerticesByEdge(newedgeinserter.first, subgraph); _subgraph_edge_index.insert( std::make_pair(subgraph_edge_by_vertices, edgeiter->second) ); LOG(4, "DEBUG: Adding edge " << sourceindex << "<->" << targetindex << " in graph to subgraph as edge " << subsourceindex << "<->" << subtargetindex << "."); } } } return subgraph; } bool compareEnergyContribution( const std::pair &a, const std::pair &b) { return a.second.energy < b.second.energy; } ActionState::ptr GraphChemicalSpaceEvaluatorAction::performCall() { /// 1. create boost::graph from graph6 string BoostGraphCreator BGC; BGC.createFromGraph6String(params.graph_string.get()); BoostGraphCreator::UndirectedGraph graph = BGC.get(); Extractors::index_map_t index_map = boost::get(boost::vertex_index, graph); 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 elementnumbers; elementnumbers.reserve(params.graph_elements.get().size()); size_t i = 0; for (std::vector::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(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_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 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."); Extractors::set_of_nodes_t set_of_nodes; for (size_t nodeid = 0; nodeid < nodes_in_graph; ++nodeid) { const size_t &rootindex = nodeid; LOG(2, "DEBUG: Current root index is " << rootindex); /// 5a. from node in graph (with this index) perform BFS till n-1 (#nodes in binding model) std::vector distances(boost::num_vertices(graph)); boost::breadth_first_search( graph, boost::vertex(rootindex, graph), boost::visitor(Extractors::record_distance(&distances[0]))); LOG(3, "DEBUG: BFS discovered the following distances " << distances); /// 5b. and store all node in map with distance to root as key Extractors::nodes_per_level_t nodes_per_level; for (size_t i=0;isecond.energy << " Ht till " << highest_contribution_graph->second.energy << " Ht, picking lowest."); total_energy += lowest_contribution_graph->second.energy; } } LOG(1, "The graph with degrees " << current_degrees << " has a total BOSSANOVA energy of " << total_energy); } LOG(1, "There have been " << num_admissible << " admissible degree combinations for the given graph."); return Action::success; } ActionState::ptr GraphChemicalSpaceEvaluatorAction::performUndo(ActionState::ptr _state) { return Action::success; } ActionState::ptr GraphChemicalSpaceEvaluatorAction::performRedo(ActionState::ptr _state){ return Action::success; } bool GraphChemicalSpaceEvaluatorAction::canUndo() { return true; } bool GraphChemicalSpaceEvaluatorAction::shouldUndo() { return true; } /** =========== end of function ====================== */