source: src/FunctionApproximation/Extractors.cpp@ 2124aa

/*
 * 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 "FunctionApproximation/Extractors.hpp"
#include "FunctionApproximation/FunctionArgument.hpp"

#include "Fragmentation/Homology/HomologyGraph.hpp"

using namespace boost::assign;

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;
}

typedef size_t level_t;
typedef size_t node_t;
typedef std::multimap< level_t, node_t > nodes_per_level_t;
typedef std::set<node_t> nodes_t;
typedef std::set<nodes_t> set_of_nodes_t;
typedef boost::property_map < boost::adjacency_list <>, boost::vertex_index_t >::type index_map_t;

typedef boost::bimap<
    boost::bimaps::set_of< size_t >,
    boost::bimaps::multiset_of< Extractors::ParticleType_t >
> type_index_lookup_t;

typedef std::set<node_t> set_type;
typedef std::set<set_type> powerset_type;

typedef boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS,
    boost::no_property, boost::no_property > UndirectedGraph;
typedef boost::subgraph< UndirectedGraph > UndirectedSubgraph;

typedef std::map< node_t, std::pair<Extractors::ParticleType_t, size_t> > node_FragmentNode_map_t;

typedef std::map< argument_t::indices_t, size_t> argument_placement_map_t;

typedef std::map<size_t, size_t> argindex_to_nodeindex_t;

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
 */
static 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;
}

static 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);
}

/**
 * I have no idea why this is so complicated with BGL ...
 *
 * This is taken from the book "The Boost Graph Library: User Guide and Reference Manual, Portable Documents",
 * chapter "Basic Graph Algorithms", example on calculating the bacon number.
 */
template <typename DistanceMap>
class distance_recorder : public boost::default_bfs_visitor
{
public:
  distance_recorder(DistanceMap dist) : d(dist) {}

  template <typename Edge, typename Graph>
  void tree_edge(Edge e, const Graph &g) const {
    typename boost::graph_traits<Graph>::vertex_descriptor u = source(e,g), v = target(e,g);
    d[v] = d[u] + 1;
  }

private:
  DistanceMap d;
};

template <typename DistanceMap>
distance_recorder<DistanceMap> record_distance(DistanceMap d)
{
  return distance_recorder<DistanceMap>(d);
}

FunctionModel::list_of_arguments_t Extractors::reorderArgumentsByParticleTypes(
    const FunctionModel::list_of_arguments_t &listargs,
    const HomologyGraph &_graph,
    const ParticleTypes_t &_types,
    const HomologyGraph &_bindingmodel
    )
{
  FunctionModel::list_of_arguments_t returnargs;
  for (FunctionModel::list_of_arguments_t::const_iterator listiter = listargs.begin();
      listiter != listargs.end(); ++listiter) {
    const FunctionModel::arguments_t &args = *listiter;

    // deal with the case when there are no distances (ConstantPotential)
    if (_bindingmodel.getNodes().size() < 2) {
      LOG(3, "DEBUG: Potential requires only zero or no particle types, needs no distances.");
      continue;
    }
    if (_bindingmodel.getEdges().empty()) {
      LOG(3, "DEBUG: Potential represents non-bonded interactions, gets all distances.");
      returnargs.push_back(args);
      continue;
    }

    /// 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 FragmentNode &firstnode = _bindingmodel.getNodes().begin()->first;
    const size_t &firsttimes = _bindingmodel.getNodes().begin()->second;
    const size_t &firsttype = firstnode.getAtomicNumber(); // TODO: convert to ParticleNumber_t ?

    /// 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);
#ifndef NDEBUG
    const size_t checktimes = std::distance(range.first, range.second);
    ASSERT( firsttimes == checktimes,
        "Extractors::reorderArgumentsByParticleTypes() - unequal amounts of same first type "
        +toString(firsttype)+" in HomologyGraph and in type lookup "+toString(checktimes));
#endif

    /// 4. go over all candidate nodes (gives index)
    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(
          N-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) {
        LOG(2, "ACCEPT: " << nodes_graph << " is identical to " << _bindingmodel);
        /// 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 store for each node its type and the number of connected
        /// edges, which we extract from the argumentbunch.
        node_FragmentNode_map_t node_FragmentNode_map = fillNodeFragmentMap(argumentbunch);

        /// Then, we step through the nodes of the bindingmodel
        const HomologyGraph::nodes_t bindingmodel_nodes = _bindingmodel.getNodes();
        argindex_to_nodeindex_t argindex_to_nodeindex;
        size_t nodeindex = 0;
        for (HomologyGraph::nodes_t::const_iterator nodeiter = bindingmodel_nodes.begin();
            nodeiter != bindingmodel_nodes.end(); ++nodeiter) {
          const FragmentNode &node = nodeiter->first;
          LOG(3, "DEBUG: Binding model's node #" << node << ".");
          // might have to pick several
          for (size_t pick_no = 0; pick_no < nodeiter->second; ++pick_no) {
            /// ... 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 as choice #"
                    << pick_no << ".");
                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);
          }
        }
        /// This gives then the desired mapping from indices in the arguments to
        /// the index in the order of the HomologyGraph's nodes

        /// 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
        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 known placement map
        FunctionModel::arguments_t sortedargs(argumentbunch);
        for (FunctionModel::arguments_t::const_iterator argiter = argumentbunch.begin();
            argiter != argumentbunch.end(); ++argiter) {
          const argument_t &arg = *argiter;
          const argument_placement_map_t::const_iterator indexiter =
              argument_placement_map.find( translateIndices(argindex_to_nodeindex, arg.indices) );
          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 " << translateIndices(argindex_to_nodeindex, arg.indices));
        }
        LOG(2, "DEBUG: Sorted arguments are " << sortedargs << ".");

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

  return returnargs;
}

FunctionModel::list_of_arguments_t Extractors::filterArgumentsByParticleTypes(
    const FunctionModel::arguments_t &args,
    const HomologyGraph &_graph,
    const ParticleTypes_t &_types,
    const HomologyGraph &_bindingmodel
    )
{
  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, boost::ref(args) > > ArgsLookup_t;
//  ArgsLookup_t ArgsLookup;

  // 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::arguments_t allargs;
  allargs.reserve(args.size());
  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)) {
          allargs.push_back( *iter );
          iter = availableList.erase(iter);
          LOG(4, "DEBUG: Accepted argument.");
        } else if ((iter->types.first == *seconditer)
              && (iter->types.second == *firstiter)) {
          allargs.push_back( *iter );
          iter = availableList.erase(iter);
          LOG(4, "DEBUG: Accepted (flipped) argument.");
        } else {
          ++iter;
          LOG(4, "DEBUG: Rejected argument.");
        }
      }
    }
  }
  LOG(2, "DEBUG: Final list of args is " << allargs << ".");

  // first, we bring together tuples of distances that belong together
  FunctionModel::list_of_arguments_t singlelist_allargs;
  singlelist_allargs.push_back(allargs);
  FunctionModel::list_of_arguments_t returnargs =
      reorderArgumentsByParticleTypes(singlelist_allargs, _graph, _types, _bindingmodel);
//  ASSERT( sortedargs.size() == (size_t)1,
//      "Extractors::filterArgumentsByParticleTypes() - reordering did not generate a single list.");
//  // then we split up the tuples of arguments and place each into single list
//  FunctionModel::list_of_arguments_t returnargs;
//  FunctionModel::arguments_t::const_iterator argiter = sortedargs.begin()->begin();
//  const size_t num_types = _types.size();
//  const size_t args_per_tuple = num_types * (num_types-1) / 2;
//  while (argiter != sortedargs.begin()->end()) {
//    FunctionModel::arguments_t currenttuple(args_per_tuple);
//    const FunctionModel::arguments_t::const_iterator startiter = argiter;
//    std::advance(argiter, args_per_tuple);
//#ifndef NDEBUG
//    FunctionModel::arguments_t::const_iterator endoutiter =
//#endif
//        std::copy(startiter, argiter, currenttuple.begin());
//    ASSERT( endoutiter == currenttuple.end(),
//        "Extractors::filterArgumentsByParticleTypes() - currenttuple not initialized to right size.");
//    returnargs.push_back(currenttuple);
//  }

  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;
}