source: src/FunctionApproximation/Extractors.cpp@ d83d64

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2012 University of Bonn. All rights reserved.
 * Copyright (C)  2013 Frederik Heber. All rights reserved.
 * Please see the COPYING file or "Copyright notice" in builder.cpp for details.
 *
 *
 *   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/>.
 */

/*
 * Extractors.cpp
 *
 *  Created on: 15.10.2012
 *      Author: heber
 */

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

#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/breadth_first_search.hpp>
#include <boost/graph/subgraph.hpp>

//#include "CodePatterns/MemDebug.hpp"

#include <algorithm>
#include <iterator>
#include <set>
#include <sstream>
#include <utility>
#include <vector>
#include <boost/assign.hpp>
#include <boost/bimap.hpp>
#include <boost/bimap/set_of.hpp>
#include <boost/bimap/multiset_of.hpp>
#include <boost/bind.hpp>
#include <boost/foreach.hpp>

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/IteratorAdaptors.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/toString.hpp"

#include "LinearAlgebra/Vector.hpp"

#include "Fragmentation/Homology/HomologyGraph.hpp"
#include "FunctionApproximation/Extractors.hpp"
#include "FunctionApproximation/FunctionArgument.hpp"
#include "Potentials/BindingModel.hpp"

using namespace boost::assign;

using namespace Extractors;

FunctionModel::arguments_t
Extractors::gatherAllSymmetricDistanceArguments(
    const Fragment::positions_t& positions,
    const Fragment::atomicnumbers_t& atomicnumbers,
    const FragmentationEdges::edges_t &edges,
    const size_t globalid)
{
  FunctionModel::arguments_t result;

  // place edges in map
  typedef std::set< std::pair<size_t, size_t> > sorted_edges_t;
  sorted_edges_t sorted_edges;
  for (FragmentationEdges::edges_t::const_iterator edgeiter = edges.begin();
      edgeiter != edges.end(); ++edgeiter) {
    std::pair<sorted_edges_t::iterator, bool> inserter =
        sorted_edges.insert(
        (edgeiter->first < edgeiter->second) ?
            std::make_pair(edgeiter->first, edgeiter->second) :
            std::make_pair(edgeiter->second, edgeiter->first));
    ASSERT(inserter.second,
        "Extractors::gatherAllSymmetricDistanceArguments() - edge "
        +toString(*inserter.first)+" is already present");
  }

  // go through current configuration and gather all other distances
  Fragment::positions_t::const_iterator firstpositer = positions.begin();
  for (;firstpositer != positions.end(); ++firstpositer) {
    Fragment::positions_t::const_iterator secondpositer = firstpositer;
    for (; secondpositer != positions.end(); ++secondpositer) {
      if (firstpositer == secondpositer)
        continue;
      argument_t arg;
      const Vector firsttemp((*firstpositer)[0],(*firstpositer)[1],(*firstpositer)[2]);
      const Vector secondtemp((*secondpositer)[0],(*secondpositer)[1],(*secondpositer)[2]);
      arg.distance = firsttemp.distance(secondtemp);
      arg.types = std::make_pair(
          (int)atomicnumbers[ std::distance(positions.begin(), firstpositer) ],
          (int)atomicnumbers[ std::distance(positions.begin(), secondpositer) ]
          );
      arg.indices = std::make_pair(
          std::distance(
              positions.begin(), firstpositer),
          std::distance(
              positions.begin(), secondpositer)
          );
      arg.globalid = globalid;
      arg.bonded = sorted_edges.find(arg.indices) != sorted_edges.end();
      LOG(3, "DEBUG: Created argument " << arg << ".");
      result.push_back(arg);
    }
  }

  return result;
}

Extractors::elementcounts_t 
Extractors::_detail::getElementCounts(
    const Fragment::atomicnumbers_t elements
  )
{
  elementcounts_t elementcounts;
  for (Fragment::atomicnumbers_t::const_iterator elementiter = elements.begin();
      elementiter != elements.end(); ++elementiter) {
    // insert new element
    std::pair< elementcounts_t::iterator, bool> inserter =
        elementcounts.insert( std::make_pair( *elementiter, 1) );
    // if already present, just increase its count
    if (!inserter.second)
      ++(inserter.first->second);
  }
  return elementcounts;
}

struct ParticleTypesComparator {
  bool operator()(const argument_t::types_t &a, const argument_t::types_t &b)
  {
    if (a.first < a.second) {
      if (b.first < b.second) {
        if (a.first < b.first)
          return true;
        else if (a.first > b.first)
          return false;
        else
          return (a.second < b.second);
      } else {
        if (a.first < b.second)
          return true;
        else if (a.first > b.second)
          return false;
        else
          return (a.second < b.first);
      }
    } else {
      if (b.first < b.second) {
        if (a.second < b.first)
          return true;
        else if (a.second > b.first)
          return false;
        else
          return (a.first < b.second);
      } else {
        if (a.second < b.second)
          return true;
        else if (a.second > b.second)
          return false;
        else
          return (a.first < b.first);
      }
    }
  }
};

std::ostream& operator<<(std::ostream &out, const argument_t::types_t &a)
{
  out << "[" << a.first << "," << a.second << "]";
  return out;
}

void insertIntoNodeFragmentMap(
    node_FragmentNode_map_t &_node_FragmentNode_map,
    const size_t &_index,
    const Extractors::ParticleType_t &_type)
{
  const node_FragmentNode_map_t::iterator mapiter = _node_FragmentNode_map.find(_index);
  // check if already present
  if (mapiter != _node_FragmentNode_map.end()) {
    // assert same type and increment number of edges
    ASSERT( mapiter->second.first == _type,
        "insertIntoNodeFragmentMap() - different types "+toString(mapiter->second.first)+
        " and "+toString(_type)+" for node "+toString(_index));
    ++(mapiter->second.second);
  } else {
    // place new entry with a single present edge
    _node_FragmentNode_map.insert( std::make_pair(_index, std::make_pair(_type, 1) ));
  }
}

static node_FragmentNode_map_t fillNodeFragmentMap(
    FunctionModel::arguments_t &argumentbunch)
{
  node_FragmentNode_map_t node_FragmentNode_map;
  // place each node index with type and number of edges into map
  for (FunctionModel::arguments_t::const_iterator argiter = argumentbunch.begin();
      argiter != argumentbunch.end(); ++argiter) {
    const argument_t &arg = *argiter;
    // only consider the distances that represent a bond edge
    if (arg.bonded) {
      insertIntoNodeFragmentMap(node_FragmentNode_map, arg.indices.first, arg.types.first);
      insertIntoNodeFragmentMap(node_FragmentNode_map, arg.indices.second, arg.types.second);
    }
  }

  return node_FragmentNode_map;
}

static argument_placement_map_t fillArgumentsPlacementMap(const size_t num_args)
{
  argument_placement_map_t argument_placement_map;
  size_t placement = 0;
  for (size_t i = 0;i<num_args;++i)
    for (size_t j = i+1;j<num_args;++j)
      argument_placement_map.insert( std::make_pair( argument_t::indices_t(i,j), placement++) );
  return argument_placement_map;
}

static argument_t::indices_t translateIndices(
    const argindex_to_nodeindex_t &_argindex_to_nodeindex,
    const argument_t::indices_t &_indices
    )
{
  const argindex_to_nodeindex_t::const_iterator firstiter =
      _argindex_to_nodeindex.find(_indices.first);
  ASSERT( firstiter != _argindex_to_nodeindex.end(),
      "translateIndices() - could not find #"+toString(_indices.first)+
      " in translation map.");
  const argindex_to_nodeindex_t::const_iterator seconditer =
      _argindex_to_nodeindex.find(_indices.second);
  ASSERT(seconditer != _argindex_to_nodeindex.end(),
      "translateIndices() - could not find #"+toString(_indices.second)+
      " in translation map.");
  // mind increasing order of indices in pair
  if (firstiter->second < seconditer->second)
    return argument_t::indices_t(firstiter->second, seconditer->second);
  else
    return argument_t::indices_t(seconditer->second, firstiter->second);
}

/** Power set generator
 *
 * taken from https://rosettacode.org/wiki/Power_set#Non-recursive_version
 *
 */
static powerset_type powerset(set_type const& set)
{
  typedef set_type::const_iterator set_iter;
  typedef std::vector<set_iter> vec;

  struct local
  {
    static int dereference(set_iter v) { return *v; }
  };

  powerset_type result;

  vec elements;
  do {
    set_type tmp;
    std::transform(elements.begin(), elements.end(),
                   std::inserter(tmp, tmp.end()),
                   local::dereference);
    result.insert(tmp);
    if (!elements.empty() && ++elements.back() == set.end()) {
      elements.pop_back();
    } else {
      set_iter iter;
      if (elements.empty()) {
        iter = set.begin();
      } else {
        iter = elements.back();
        ++iter;
      }
      for (; iter != set.end(); ++iter) {
        elements.push_back(iter);
      }
    }
  } while (!elements.empty());

  return result;
}

/** Recursive function to generate all induced, connected subgraphs given a
 * graph.
 *
 * \param N number of still left to pick
 * \param depth level in \a set_of_nodes
 * \param nodes current set of nodes that are picked already
 * \param set_of_nodes resulting set of generated subgraphs' nodes
 * \param nodes_per_level level-wise frontier of connected nodes around a root node
 * \param graph graph containing the adjacency
 * \param index_map with indices per \a graph' vertex
 */
void generateAllInducedConnectedSubgraphs(
    const size_t N,
    const level_t level,
    const nodes_t &nodes,
    set_of_nodes_t &set_of_nodes,
    const nodes_per_level_t &nodes_per_level,
    const UndirectedGraph &graph,
    const std::vector<size_t> &_distance,
    const index_map_t &index_map)
{
  ASSERT( nodes_per_level.find(level) != nodes_per_level.end(),
      "generateAllInducedConnectedSubgraphs() - we are deeper than the graph.");
  ASSERT( N < nodes_per_level.size(),
      "generateAllInducedConnectedSubgraphs() - we are looking for subgraphs larger than the graph.");
  if (N > 0) {
    LOG(3, "DEBUG: At level " << level << " current nodes is " << nodes << ", need to find " << N << " more.");
    // get next level's set and constrain to nodes connected to this set
    nodes_t validnodes;
    std::pair< nodes_per_level_t::const_iterator, nodes_per_level_t::const_iterator> range =
        nodes_per_level.equal_range(level);
    for (nodes_per_level_t::const_iterator rangeiter = range.first;
        rangeiter != range.second; ++rangeiter) {
      LOG(4, "DEBUG: Checking edges further away from node #" << rangeiter->second);
      // get all edges connected to this node further away
      UndirectedGraph::in_edge_iterator i, end;
      boost::tie(i, end) = boost::in_edges(boost::vertex(rangeiter->second, graph), graph);
      for (;i != end; ++i) {
        // check each edge whether it's in nodes
        const node_t sourceindex = boost::get(index_map, boost::source(*i, graph));
        const node_t targetindex = boost::get(index_map, boost::target(*i, graph));
        const size_t &source_distance = _distance[sourceindex];
        const size_t &target_distance = _distance[targetindex];
        // edge is going deeper into graph
        if (((source_distance == level) && (target_distance == (level+1)))
            || ((source_distance == (level+1)) && (target_distance == level))) {
          LOG(5, "DEBUG: Candidate edge on level " << level << " is from " << sourceindex
              << " to " << targetindex << ".");
          // pick right index and check for not present in list yet
          if (sourceindex == rangeiter->second) {
            if (nodes.count(targetindex) == 0) {
              validnodes.insert(targetindex);
              LOG(4, "DEBUG: Inserting node #" << targetindex << " into valid nodes.");
            }
          } else if (targetindex == rangeiter->second) {
            if (nodes.count(sourceindex) == 0) {
              validnodes.insert(sourceindex);
              LOG(4, "DEBUG: Inserting node #" << sourceindex << " into valid nodes.");
            }
          } else {
            ASSERT(0,
                "generateAllInducedConnectedSubgraphs() - neither source #"+toString(sourceindex)+
                " nor target #"+toString(targetindex)+" is equal to #"+toString(rangeiter->second));
          }
        }
      }
    }
    // skip this if we cannot go deeper into the graph from here
    if (validnodes.empty()) {
      LOG(3, "DEBUG: We did not find any more nodes to step on from " << nodes << ".");
      return;
    }

    // go through power set
    const powerset_type test_powerset = powerset(validnodes);
    for (powerset_type::const_iterator iter = test_powerset.begin();
         iter != test_powerset.end();
         ++iter) {
      // count set bits (#elements in *iter), must be between 1 and N
      const size_t num_set_bits = iter->size();
      if ((num_set_bits > 0) && (num_set_bits <= N)) {
        // add set to nodes
        LOG(3, "DEBUG: With present " << nodes << " the current set is " << *iter << " of "
            << validnodes << ".");

        // copy the nodes before insertion
        nodes_t filled_nodes(nodes.begin(), nodes.end());
        filled_nodes.insert(iter->begin(), iter->end());

        // and recurse
        generateAllInducedConnectedSubgraphs(
            N-num_set_bits, level+1, filled_nodes, set_of_nodes, nodes_per_level, graph, _distance, index_map);
      }
    }
  } else {
    // N==0: we have a winner
    std::pair<set_of_nodes_t::iterator, bool> inserter =
        set_of_nodes.insert( nodes );
    if (!inserter.second)
      LOG(2, "DEBUG: subgraph " << nodes << " is already contained in set_of_nodes.");
    else
      LOG(2, "DEBUG: subgraph " << nodes << " inserted into set_of_nodes.");
  }
}

static Extractors::ParticleType_t getParticleTypeToNode(
    const type_index_lookup_t &type_index_lookup,
    const size_t nodeindex)
{
  const type_index_lookup_t::left_const_iterator typeiter = type_index_lookup.left.find(nodeindex);
  ASSERT( typeiter != type_index_lookup.left.end(),
      "getParticleTypeToNode() - could not find type to node #"+toString(nodeindex));
  return typeiter->second;
}

HomologyGraph createHomologyGraphFromNodes(
    const nodes_t &nodes,
    const type_index_lookup_t &type_index_lookup,
    const UndirectedGraph &graph,
    const index_map_t &index_map
    )
{
  HomologyGraph::nodes_t graph_nodes;
  HomologyGraph::edges_t graph_edges;
  {
    typedef std::set< std::pair<node_t, node_t> > graph_edges_t;
    graph_edges_t edge_set;
    std::pair<HomologyGraph::nodes_t::iterator, bool > inserter;
    for (nodes_t::const_iterator nodeiter = nodes.begin();
        nodeiter != nodes.end(); ++nodeiter) {
      const Extractors::ParticleType_t &nodetype = getParticleTypeToNode(type_index_lookup, *nodeiter);

      // count edges in constrained set for this particular node
      size_t num_edges = 0;
      UndirectedGraph::in_edge_iterator i, end;
      for (boost::tie(i, end) = boost::in_edges(boost::vertex(*nodeiter, graph), graph);
          i != end; ++i) {
        const node_t sourceindex = boost::get(index_map, boost::source(*i, graph));
        const node_t targetindex = boost::get(index_map, boost::target(*i, graph));
        // check each edge whether it's in nodes
        if ((nodes.count(sourceindex) != 0) && (nodes.count(targetindex) != 0)) {
          // increase edge count of node
          ++num_edges;
          // we first store edges in a set to ensure their uniqueness (as we encounter
          // each edge two times and we don't know if source and target will be
          // different the second time encountered)
          if (sourceindex < targetindex)
            edge_set.insert( std::make_pair(sourceindex, targetindex) );
          else
            edge_set.insert( std::make_pair(targetindex, sourceindex) );
        }
      }
      LOG(4, "DEBUG: Node #" << *nodeiter << " has " << num_edges << " edges.");

      // add FragmentNode
      inserter = graph_nodes.insert( std::make_pair(FragmentNode(nodetype, num_edges), 1) );
      if (!inserter.second)
        ++(inserter.first->second);
    }

    // add edges
    for (graph_edges_t::const_iterator edgeiter = edge_set.begin();
        edgeiter != edge_set.end(); ++edgeiter) {
      const Extractors::ParticleType_t sourcetype =
          getParticleTypeToNode(type_index_lookup, edgeiter->first);
      const Extractors::ParticleType_t targettype =
          getParticleTypeToNode(type_index_lookup, edgeiter->second);
      // FragmentEdge takes care of proper sorting
      FragmentEdge edge(sourcetype, targettype);
      LOG(4, "DEBUG: Adding fragment edge " << edge);
      std::pair<HomologyGraph::edges_t::iterator, bool > inserter;
      inserter = graph_edges.insert( std::make_pair( edge, 1) );
      if (!inserter.second)
        ++(inserter.first->second);
    }
  }

  return HomologyGraph(graph_nodes, graph_edges);
}

FunctionModel::list_of_arguments_t Extractors::filterArgumentsByBindingModel(
    const FunctionModel::arguments_t &args,
    const HomologyGraph &_graph,
    const ParticleTypes_t &_types,
    const BindingModel &_bindingmodel
    )
{
  FunctionModel::list_of_arguments_t returnargs;

  // deal with the case when there are no distances (ConstantPotential)
  if (_bindingmodel.getNodes().size() < 2) {
    LOG(3, "DEBUG: Potential requires only one or no particle types, needs no distances.");
    return returnargs;
  }
  if (_bindingmodel.getGraph().getEdges().empty()) {
    LOG(3, "DEBUG: Potential represents non-bonded interactions, gets all distances.");
    // TODO: Here we need to constrain to all distances matching the types?
    returnargs.push_back(args);
    return returnargs;
  }

  /// 0. place all nodes in a lookup map for their type
  type_index_lookup_t type_index_lookup;
  for(FunctionModel::arguments_t::const_iterator iter = args.begin();
      iter != args.end(); ++iter) {
    if (type_index_lookup.left.find(iter->indices.first) == type_index_lookup.left.end())
      type_index_lookup.left.insert( std::make_pair(iter->indices.first, iter->types.first) );
    else {
      ASSERT(type_index_lookup.left.at(iter->indices.first) == iter->types.first,
          "Extractors::reorderArgumentsByParticleTypes() - entry " +toString(iter->indices.first)+
          " is already present with different type "+toString(iter->types.first)
          +" than "+toString(type_index_lookup.left.at(iter->indices.first)));
    }
    if (type_index_lookup.left.find(iter->indices.second) == type_index_lookup.left.end())
      type_index_lookup.left.insert( std::make_pair(iter->indices.second, iter->types.second) );
    else {
      ASSERT(type_index_lookup.left.at(iter->indices.second) == iter->types.second,
          "Extractors::reorderArgumentsByParticleTypes() - entry " +toString(iter->indices.second)+
          " is already present with different type "+toString(iter->types.second)
          +" than "+toString(type_index_lookup.left.at(iter->indices.second)));
    }
  }
  if (DoLog(3)) {
    std::stringstream output;
    for (type_index_lookup_t::left_const_iterator indexiter = type_index_lookup.left.begin();
        indexiter != type_index_lookup.left.end(); ++indexiter) {
      output << " [" << indexiter->first << "," << indexiter->second << "]";
    }
    LOG(3, "DEBUG: index to type map:" << output.str());
  }
  if (DoLog(3)) {
    std::stringstream output;
    for (type_index_lookup_t::right_const_iterator typeiter = type_index_lookup.right.begin();
        typeiter != type_index_lookup.right.end(); ++typeiter) {
      output << " [" << typeiter->first << "," << typeiter->second << "]";
    }
    LOG(3, "DEBUG: type to index map " << output.str());
  }

  /// 1. construct boost::graph from arguments_t (iter)
  const size_t N = type_index_lookup.left.size();
  UndirectedGraph fragmentgraph(N);
  for(FunctionModel::arguments_t::const_iterator iter = args.begin();
      iter != args.end(); ++iter) {
    if (iter->bonded)
      boost::add_edge(iter->indices.first, iter->indices.second, fragmentgraph);
  }
  index_map_t index_map = boost::get(boost::vertex_index, fragmentgraph);
  LOG(2, "DEBUG: We have " << boost::num_vertices(fragmentgraph) << " nodes in the fragment graph.");

  /// 2. grab first node of the binding model (gives a type)
  const BindingModel::vector_nodes_t::const_iterator firstiter = _bindingmodel.getNodes().begin();
  const FragmentNode &firstnode = *firstiter;
  const Extractors::ParticleType_t &firsttype = firstnode.getAtomicNumber();

  /// 3. grab all candidate nodes contained in arguments_t
  std::pair<
      type_index_lookup_t::right_const_iterator,
      type_index_lookup_t::right_const_iterator> range = type_index_lookup.right.equal_range(firsttype);

  /// 4. go over all candidate nodes (gives index)
  const size_t nodes_in_bindingmodel = _bindingmodel.getNodes().size();
  LOG(2, "DEBUG: There are " << nodes_in_bindingmodel << " nodes in the binding model.");
  set_of_nodes_t set_of_nodes;
  for (type_index_lookup_t::right_const_iterator rangeiter = range.first;
      rangeiter != range.second; ++rangeiter) {
    const size_t &rootindex = rangeiter->second;
    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(fragmentgraph));
    boost::breadth_first_search(
        fragmentgraph,
        boost::vertex(rootindex, fragmentgraph),
        boost::visitor(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
    nodes_per_level_t nodes_per_level;
    for (size_t i=0;i<distances.size();++i) {
      nodes_per_level.insert( std::make_pair(level_t(distances[i]), node_t(i)) );
    }
    LOG(3, "DEBUG: Nodes per level are " << nodes_per_level);

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

    generateAllInducedConnectedSubgraphs(
        nodes_in_bindingmodel-1, level, nodes, set_of_nodes, nodes_per_level, fragmentgraph, distances, index_map);
    LOG(2, "DEBUG: We have found " << set_of_nodes.size() << " unique induced, connected subgraphs.");
  }

  /// 8. go through each induced, connected subgraph
  for (set_of_nodes_t::const_iterator nodesiter = set_of_nodes.begin();
      nodesiter != set_of_nodes.end(); ++nodesiter) {
    /// 9. convert its nodes into a HomologyGraph using the subgraph and arguments_t (types for index)
    const nodes_t &nodes = *nodesiter;
    const HomologyGraph nodes_graph =
        createHomologyGraphFromNodes(nodes, type_index_lookup, fragmentgraph, index_map);

    /// 10. compare each converted HomologyGraph with _bindingmodel: Accept or Reject
    if (nodes_graph == _bindingmodel.getGraph()) {
      LOG(2, "ACCEPT: " << nodes_graph << " is identical to " << _bindingmodel.getGraph());
      /// 11. for each accepted keyset, pick _all_ symmetric distances from arguments_t
      FunctionModel::arguments_t argumentbunch;
      argumentbunch.reserve(args.size());
      for(FunctionModel::arguments_t::const_iterator iter = args.begin();
          iter != args.end(); ++iter) {
        if ((nodes.count(iter->indices.first) != 0) && (nodes.count(iter->indices.second) != 0)) {
          argumentbunch.push_back(*iter);
        }
      }
      const size_t num_symmetric_distances = nodes.size()*(nodes.size()-1)/2;
      ASSERT( argumentbunch.size() == num_symmetric_distances,
          "Extractors::reorderArgumentsByParticleTypes() - only "+toString(argumentbunch.size())+
          " found instead of "+toString(num_symmetric_distances));
      LOG(3, "DEBUG: Final bunch of arguments is " << argumentbunch << ".");

      /**
       * We still need to bring the arguments in the correct order for the potential.
       *
       * The potential gives the desired order via the nodes in the HomologyGraph! Their
       * sequence matches with the user-defined particle types and because of the properties
       * of the homology graph (FragmentNode stores type and number of edges) it also
       * matches with the desired bonding graph.
       */

      /// So, we need to extract from the argumentbunch the information contained in each
      /// FragmentNode, namely its type and the number of connected edges
      node_FragmentNode_map_t node_FragmentNode_map = fillNodeFragmentMap(argumentbunch);

      /// Then, we step through the nodes of the bindingmodel ...
      /// ... and find a suitable mapping from indices in the arguments to
      /// the index in the order of the HomologyGraph's nodes
      const BindingModel::vector_nodes_t bindingmodel_nodes = _bindingmodel.getNodes();
      argindex_to_nodeindex_t argindex_to_nodeindex;
      size_t nodeindex = 0;
      for (BindingModel::vector_nodes_t::const_iterator nodeiter = bindingmodel_nodes.begin();
          nodeiter != bindingmodel_nodes.end(); ++nodeiter) {
        const FragmentNode &node = *nodeiter;
        LOG(3, "DEBUG: Binding model's node #" << node << ".");
        /// ... and pick for each the first (and unique) from these stored nodes.
        node_FragmentNode_map_t::iterator mapiter = node_FragmentNode_map.begin();
        for (;mapiter != node_FragmentNode_map.end(); ++mapiter) {
          if ((mapiter->second.first == node.getAtomicNumber())
              && (mapiter->second.second == node.getConnectedEdges())) {
            LOG(3, "DEBUG: #" << mapiter->first << " with type " << mapiter->second.first
                << " and " << mapiter->second.second << " connected edges matches.");
            break;
          }
        }
        ASSERT( mapiter != node_FragmentNode_map.end(),
            "Extractors::reorderArgumentsByParticleTypes() - could not find a suitable node for #"+
            toString(mapiter->first)+" with type "+toString(mapiter->second.first)+" and "+
            toString(mapiter->second.second)+" connected edges");
        std::pair<argindex_to_nodeindex_t::iterator, bool> inserter =
            argindex_to_nodeindex.insert( std::make_pair(mapiter->first, nodeindex++) );
        ASSERT( inserter.second,
            "Extractors::reorderArgumentsByParticleTypes() - node #"+toString(mapiter->first)+
            " is already present?");
        // remove to ensure uniqueness of choice
        node_FragmentNode_map.erase(mapiter);
      }
      LOG(4, "DEBUG: argument's indices to node index map is " << argindex_to_nodeindex);
      // i.e. this is not the arg's index in argumentbunch, but the index of the position
      // contained in the argument_t

      /// Finally, we only need to bring the arguments in the typical order:
      /// 01 02 03 04 ... 0n, 12 13 14 ... 1n, 23 24 25 ... 2n, ...,  (n-1)n
      /// These ordering we store in a map for each argument's indices
      const size_t num_args = argindex_to_nodeindex.size();
      argument_placement_map_t argument_placement_map = fillArgumentsPlacementMap(num_args);
      LOG(4, "DEBUG: Placement map is " << argument_placement_map);
      ASSERT( argument_placement_map.size() == argumentbunch.size(),
          "Extractors::reorderArgumentsByParticleTypes() - placement map has size "+
          toString(argument_placement_map.size())+" and we expected "+toString(argumentbunch.size()));

      // and finally resort the arguments with the constructed placement map
      FunctionModel::arguments_t sortedargs(argumentbunch.size());
      for (FunctionModel::arguments_t::const_iterator argiter = argumentbunch.begin();
          argiter != argumentbunch.end(); ++argiter) {
        const argument_t &arg = *argiter;
        const argument_t::indices_t translatedIndices =
            translateIndices(argindex_to_nodeindex, arg.indices);
        const argument_placement_map_t::const_iterator indexiter =
            argument_placement_map.find( translatedIndices );
        ASSERT( indexiter != argument_placement_map.end(),
            "Extractors::reorderArgumentsByParticleTypes() - could not find place for edge "+
            toString(arg.indices));
        sortedargs[indexiter->second] = arg;
        LOG(3, "DEBUG: Placed argument " << arg << " at #" << indexiter->second
            << " as translated indices are " << translatedIndices);
      }
      LOG(2, "DEBUG: Sorted arguments are " << sortedargs << ".");

      returnargs.push_back(sortedargs);
    } else {
      LOG(2, "REJECT: " << nodes_graph << " is not identical to " << _bindingmodel.getGraph());
    }
  }

  return returnargs;
}

FunctionModel::list_of_arguments_t Extractors::filterArgumentsByParticleTypes(
    const FunctionModel::arguments_t &args,
    const HomologyGraph &_graph,
    const ParticleTypes_t &_types,
    const BindingModel &_bindingmodel
    )
{
  // list allows for quicker removal than vector
  typedef std::list< argument_t > ListArguments_t;
  ListArguments_t availableList(args.begin(), args.end());
  LOG(2, "DEBUG: Initial list of args is " << args << ".");

  // TODO: fill a lookup map such that we don't have O(M^3) scaling, if M is number
  // of types (and we always must have M(M-1)/2 args) but O(M^2 log(M)). However, as
  // M is very small (<=3), this is not necessary fruitful now.
//  typedef ParticleTypes_t firsttype;
//  typedef ParticleTypes_t secondtype;
//  typedef std::map< firsttype, std::map< secondtype, FunctionModel::arguments_t > > ArgsLookup_t;
//  ArgsLookup_t ArgsLookup;

  ASSERT( _types.size() <= 2,
      "Extractors::filterArgumentsByParticleTypes() - this only filters and does not reorder."
      +std::string(" Hence, it is not useful for multiple arguments per model."));

  // basically, we have two choose any two pairs out of types but only those
  // where the first is less than the latter. Hence, we start the second
  // iterator at the current position of the first one and skip the equal case.
  FunctionModel::list_of_arguments_t returnargs;
  for (ParticleTypes_t::const_iterator firstiter = _types.begin();
      firstiter != _types.end();
      ++firstiter) {
    for (ParticleTypes_t::const_iterator seconditer = firstiter;
        seconditer != _types.end();
        ++seconditer) {
      if (seconditer == firstiter)
        continue;
      LOG(3, "DEBUG: Looking for (" << *firstiter << "," << *seconditer << ") in all args.");

      // search the right one in _args (we might allow switching places of
      // firstiter and seconditer, as distance is symmetric).
      ListArguments_t::iterator iter = availableList.begin();
      while (iter != availableList.end()) {
        LOG(4, "DEBUG: Current args is " << *iter << ".");
        if ((iter->types.first == *firstiter)
              && (iter->types.second == *seconditer)) {
          returnargs.push_back( FunctionModel::arguments_t(1, *iter) );
          iter = availableList.erase(iter);
          LOG(4, "DEBUG: Accepted argument.");
        } else if ((iter->types.first == *seconditer)
              && (iter->types.second == *firstiter)) {
          returnargs.push_back( FunctionModel::arguments_t(1, *iter) );
          iter = availableList.erase(iter);
          LOG(4, "DEBUG: Accepted (flipped) argument.");
        } else {
          ++iter;
          LOG(4, "DEBUG: Rejected argument.");
        }
      }
    }
  }

  LOG(2, "DEBUG: We have generated " << returnargs.size() << " tuples of distances.");

  return returnargs;
}


FunctionModel::arguments_t Extractors::combineArguments(
    const FunctionModel::arguments_t &firstargs,
    const FunctionModel::arguments_t &secondargs)
{
  FunctionModel::arguments_t args = concatenateArguments(firstargs, secondargs);
  std::sort(args.begin(), args.end(),
      boost::bind(&argument_t::operator<, _1, _2));
  FunctionModel::arguments_t::iterator iter =
      std::unique(args.begin(), args.end(), 
          boost::bind(&argument_t::operator==, _1, _2));
  args.erase(iter, args.end());
  return args;
}

FunctionModel::arguments_t Extractors::concatenateArguments(
    const FunctionModel::arguments_t &firstargs,
    const FunctionModel::arguments_t &secondargs)
{
  FunctionModel::arguments_t args(firstargs);
  args.insert(args.end(), secondargs.begin(), secondargs.end());
  return args;
}

FunctionModel::list_of_arguments_t Extractors::concatenateListOfArguments(
    const FunctionModel::list_of_arguments_t &firstlistargs,
    const FunctionModel::list_of_arguments_t &secondlistargs)
{
  FunctionModel::list_of_arguments_t listargs(firstlistargs);
  listargs.insert(listargs.end(), secondlistargs.begin(), secondlistargs.end());
  return listargs;
}