source: src/Actions/GraphAction/ChemicalSpaceEvaluatorAction.cpp@ 1be100

/*
 * 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 <http://www.gnu.org/licenses/>.
 */

/*
 * ChemicalSpaceEvaluator.cpp
 *
 *  Created on: Sep 26, 2017
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

//#include "CodePatterns/MemDebug.hpp"

#include <iostream>
#include <list>
#include <map>
#include <string>
#include <vector>

#include <boost/graph/copy.hpp>

#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 <boost/graph/adjacency_list.hpp>

using namespace MoleCuilder;

// and construct the stuff
#include "ChemicalSpaceEvaluatorAction.def"
#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 = 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 <class graph_type>
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 &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);
      }
    }
  }

  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<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 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_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);
      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 <typename graph_t>
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<Extractors::node_t, BoostGraphCreator::Vertex> 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<BoostGraphCreator::Edge, bool> 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;
}


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<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(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.");
  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<size_t> 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;i<distances.size();++i) {
      nodes_per_level.insert( std::make_pair(Extractors::level_t(distances[i]), Extractors::node_t(i)) );
    }
    LOG(3, "DEBUG: Nodes per level are " << nodes_per_level);

    /// 5c. construct all possible induced connected subgraphs with this map (see fragmentation)
    Extractors::nodes_t nodes;
    // rootindex is always contained
    nodes.insert( rootindex );
    Extractors::level_t level = 0;

    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.");
  }

  // put all nodes in a list
  Extractors::nodes_t nodes;
  for (size_t nodeid = 0; nodeid < nodes_in_graph; ++nodeid) {
    nodes.insert(nodeid);
  }

  // 6. for every combination saturate fragments for lookup and energy contribution summation
  size_t num_admissible = 0;
  const HomologyContainer &container = World::getInstance().getHomologies();
  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;
    LOG(2, "DEBUG: Current degree combination is " << current_degrees);

    // saturate subgraph: add to nodes, add to type_index_lookup, index_map
      // saturate any leftover valence with hydrogen nodes
    SaturatedGraph fullgraph =
        saturateGraph(graph, 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) {
      LOG(2, "DEBUG: The degree combination " << current_degrees << " can NOT be fulfilled.");
      continue;
    } else {
      ++num_admissible;
      LOG(2, "DEBUG: The degree combination is admissable.");
    }


    // 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: Fragment nodes for graph are " << current_nodes);

      // create subgraph
      edge_index_t subgraph_edge_index;
      Extractors::UndirectedGraph newgraph = extractSubgraph(
          graph, current_nodes, edge_index, subgraph_edge_index);

      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 in fragment has index " << vindex);
        }
      }

      // saturate subgraph: add to nodes, add to type_index_lookup, index_map
        // saturate any leftover valence with hydrogen nodes
      SaturatedGraph fragmentgraph =
          saturateGraph(newgraph, new_name_map,
              type_index_lookup, type_valence_map, subgraph_edge_index,
              edges, current_degrees);

      // don't add if one of the node's valence cannot be fulfilled
      ASSERT(fragmentgraph.ValencesAllFulfilled,
          "GraphChemicalSpaceEvaluatorAction::performCall() - encountered subgraph whose valences cannot be fulfilled.");

      const Extractors::index_map_t subgraph_index_map =
          boost::get(boost::vertex_index, fragmentgraph.saturated_graph);
      // const BoostGraphCreator::name_map_t subgraph_name_map =
      //    boost::get(boost::vertex_name, fragmentgraph.saturated_graph);

      Extractors::nodes_t saturated_nodes;
      BoostGraphCreator::vertex_iter viter, vend;
      boost::tie(viter, vend) = boost::vertices(fragmentgraph.saturated_graph);
      for (; viter != vend; ++viter) {
        const Extractors::node_t vindex = boost::get(subgraph_index_map, *viter);
        LOG(4, "DEBUG: Adding fragment node index " << vindex << " of type "
            << Extractors::getParticleTypeToNode(fragmentgraph.type_index_lookup, vindex)
            << " to saturated fragment 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,
              fragmentgraph.type_index_lookup,
              fragmentgraph.saturated_graph,
              subgraph_index_map);
      LOG(2, "DEBUG: Looking for fragment graphs " << nodes_graph << " in homology container.");

      /// 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 " << 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 ====================== */