Changeset 2124aa for src/FunctionApproximation

Timestamp:

Nov 11, 2016, 2:25:30 PM (8 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_StructOpt_integration_tests, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_ChronosMutex, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, GeometryObjects, Gui_displays_atomic_force_velocity, IndependentFragmentGrids_IntegrationTest, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, RotateToPrincipalAxisSystem_UndoRedo, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, ThirdParty_MPQC_rebuilt_buildsystem, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, Ubuntu_1604_changes, stable

Children:

c738f1

Parents:

228340

git-author:

Frederik Heber <heber@…> (10/04/16 14:05:46)

git-committer:

Frederik Heber <heber@…> (11/11/16 14:25:30)

Message:

Extractors::reorderArgumentsByParticleTypes() rewritten to scan for all induced, connected subgraphs.

• after filtering the arguments down to the required particle types, we use the adjacency graph that is passed down from creating the SaturatedFragments, stored in various FragmentationResultContainers and finally used to state whether in Extractors::gatherAllSymmetricDistances() a distance represents a bond or not. In this graph we look for all subgraphs that are homologous to the graph specified by the derived EmpiricalPotential. For each accepted subgraph we then gather again all symmetric distances and concatenate them all.
• TESTFIX: Marked all potential fitting regression tests as XFAIL for the moment. On the one hand the homology container format changed (edges), and on the other hand we are in the middle of refactoring above function.
• added test using edges for ExtractorsUnitTest.
• lib dependency fixed for ..PotentialUnitTest.
• NOTE: boost::graph is incompatible with CodePattern's MEMDEBUG. Hence, includes are placed prior to memdebug's.

Location:

src/FunctionApproximation

Files:

• Extractors.cpp (modified) (6 diffs)
• unittests/ExtractorsUnitTest.cpp (modified) (1 diff)
• unittests/ExtractorsUnitTest.hpp (modified) (2 diffs)
• unittests/Makefile.am (modified) (1 diff)

src/FunctionApproximation/Extractors.cpp

@@ -35 +35 @@
 #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>
@@ -66 +76 @@
 {
   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
@@ -89 +114 @@
           );
       arg.globalid = globalid;
-      arg.bonded = false;
+      arg.bonded = sorted_edges.find(arg.indices) != sorted_edges.end();
       LOG(3, "DEBUG: Created argument " << arg << ".");
       result.push_back(arg);
@@ -161 +186 @@
 }
 
+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,
@@ -169 +539 @@
 {
   FunctionModel::list_of_arguments_t returnargs;
-  for (FunctionModel::list_of_arguments_t::const_iterator iter = listargs.begin();
-      iter != listargs.end(); ++iter) {
-    const FunctionModel::arguments_t &args = *iter;
-    /// We  place all arguments into multimap according to particle type pair.
-    // here, we need a special comparator such that types in key pair are always
-    // properly ordered.
-    typedef std::multimap<
-        argument_t::types_t,
-        argument_t,
-        ParticleTypesComparator> TypePair_Argument_Map_t;
-    TypePair_Argument_Map_t argument_map;
+  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) {
-      argument_map.insert( std::make_pair(iter->types, *iter) );
-    }
-    LOG(4, "DEBUG: particle_type map is " << argument_map << ".");
-
-    /// Then, we create the desired unique keys
-    typedef std::vector<argument_t::types_t> UniqueTypes_t;
-    UniqueTypes_t UniqueTypes;
-    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;
-        UniqueTypes.push_back( std::make_pair(*firstiter, *seconditer) );
-      }
-    }
-    LOG(4, "DEBUG: Created unique types as keys " << UniqueTypes << ".");
-
-    /// Finally, we use the unique key list to pick corresponding arguments from the map
-    FunctionModel::arguments_t sortedargs;
-    sortedargs.reserve(args.size());
-    while (!argument_map.empty()) {
-      // note that particle_types_t may be flipped, i.e. 1,8 is equal to 8,1, but we
-      // must maintain the correct order in indices in accordance with the order
-      // in _types, i.e. 1,8,1 must match with e.g. ids 1,0,2 where 1 has type 1,
-      // 0 has type 8, and 2 has type 2.
-      // In other words: We do not want to flip/modify arguments such that they match
-      // with the specific type pair we seek but then this comes at the price that we
-      // have flip indices when the types in a pair are flipped.
-
-      typedef std::vector<size_t> indices_t;
-      //!> here, we gather the indices as we discover them
-      indices_t indices;
-      indices.resize(_types.size(), (size_t)-1);
-
-      // these are two iterators that create index pairs in the same way as we have
-      // created type pairs. If a -1 is still present in indices, then the index is
-      // still arbitrary but is then set by the next found index
-      indices_t::iterator firstindex = indices.begin();
-      indices_t::iterator secondindex = firstindex+1;
-
-      //!> here, we gather the current bunch of arguments as we find them
-      FunctionModel::arguments_t argumentbunch;
-      argumentbunch.reserve(UniqueTypes.size());
-
-      for (UniqueTypes_t::const_iterator typeiter = UniqueTypes.begin();
-          typeiter != UniqueTypes.end(); ++typeiter) {
-        // have all arguments to same type pair as list within the found range
-        std::pair<
-            TypePair_Argument_Map_t::iterator,
-            TypePair_Argument_Map_t::iterator> range_t =
-                argument_map.equal_range(*typeiter);
-        LOG(4, "DEBUG: Set of arguments to current key [" << typeiter->first << ","
-            << typeiter->second << "] is " << std::list<argument_t>(
-                MapValueIterator<TypePair_Argument_Map_t::iterator>(range_t.first),
-                MapValueIterator<TypePair_Argument_Map_t::iterator>(range_t.second)
-                ) << ".");
-        // the first key is always easy and is pivot which the rest has to be associated to
-        if (typeiter == UniqueTypes.begin()) {
-          const argument_t & arg = range_t.first->second;
-          if ((typeiter->first == arg.types.first) && (typeiter->second == arg.types.second)) {
-            // store in correct order
-            *firstindex = arg.indices.first;
-            *secondindex = arg.indices.second;
-          } else {
-            // store in flipped order
-            *firstindex = arg.indices.second;
-            *secondindex = arg.indices.first;
+      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);
           }
-          argumentbunch.push_back(arg);
-          argument_map.erase(range_t.first);
-          LOG(4, "DEBUG: Gathered first argument " << arg << ".");
-        } else {
-          // go through the range and pick the first argument matching the index constraints
-          for (TypePair_Argument_Map_t::iterator argiter = range_t.first;
-              argiter != range_t.second; ++argiter) {
-            // seconditer may be -1 still
-            const argument_t &arg = argiter->second;
-            if (arg.indices.first == *firstindex) {
-              if ((arg.indices.second == *secondindex) || (*secondindex == (size_t)-1)) {
-                if (*secondindex == (size_t)-1)
-                  *secondindex = arg.indices.second;
-                argumentbunch.push_back(arg);
-                argument_map.erase(argiter);
-                LOG(4, "DEBUG: Gathered another argument " << arg << ".");
-                break;
-              }
-            } else if ((arg.indices.first == *secondindex) || (*secondindex == (size_t)-1)) {
-              if (arg.indices.second == *firstindex) {
-                if (*secondindex == (size_t)-1)
-                  *secondindex = arg.indices.first;
-                argumentbunch.push_back(arg);
-                argument_map.erase(argiter);
-                LOG(4, "DEBUG: Gathered another (flipped) argument " << arg << ".");
+        }
+        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);
           }
         }
-        // move along in indices and check bounds
-        ++secondindex;
-        if (secondindex == indices.end()) {
-          ++firstindex;
-          if (firstindex != indices.end()-1)
-            secondindex = firstindex+1;
+        /// 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));
         }
-      }
-      ASSERT( (firstindex == indices.end()-1) && (secondindex == indices.end()),
-          "Extractors::reorderArgumentsByParticleTypes() - we have not gathered enough arguments.");
-      ASSERT( argumentbunch.size() == UniqueTypes.size(),
-          "Extractors::reorderArgumentsByParticleTypes() - we have not gathered enough arguments.");
-      // place bunch of arguments in return args
-      LOG(3, "DEBUG: Given types " << _types << " and found indices " << indices << ".");
-      LOG(3, "DEBUG: Final bunch of arguments is " << argumentbunch << ".");
-      sortedargs.insert(sortedargs.end(), argumentbunch.begin(), argumentbunch.end());
-    }
-    returnargs.push_back(sortedargs);
+        LOG(2, "DEBUG: Sorted arguments are " << sortedargs << ".");
+
+        returnargs.push_back(sortedargs);
+      } else {
+        LOG(2, "REJECT: " << nodes_graph << " is not identical to " << _bindingmodel);
+      }
+    }
   }
 
@@ -368 +825 @@
   FunctionModel::list_of_arguments_t singlelist_allargs;
   singlelist_allargs.push_back(allargs);
-  FunctionModel::list_of_arguments_t sortedargs =
+  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);
-  }
+//  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.");

src/FunctionApproximation/unittests/ExtractorsUnitTest.cpp

@@ -98 +98 @@
     CPPUNIT_ASSERT( (args[i].distance >= 0) && (args[i].distance <= 4));
 }
+
+/** UnitTest for gatherAllSymmetricDistances() with edges
+ */
+void ExtractorsTest::gatherAllSymmetricDistances_edgesTest()
+{
+  // create positions
+  Fragment::positions_t positions;
+  Fragment::position_t pos(3, 0.);
+  for (double i = 0; i < 5; i+=1.) {
+    pos[0] = i;
+    positions.push_back(pos);
+  }
+  FragmentationEdges::edges_t edges;
+  for (size_t i=1;i<5;++i)
+    edges += FragmentationEdges::edge_t(i-1, i);
+
+  // create charges
+  Fragment::atomicnumbers_t atomicnumbers;
+  atomicnumbers += 6, 6, 1, 1, 1;
+
+  // create distances
+  FunctionModel::arguments_t args =
+      Extractors::gatherAllSymmetricDistances(positions, atomicnumbers, edges, 0);
+  CPPUNIT_ASSERT_EQUAL( (size_t)(5*4/2), args.size() );
+
+  // check created args
+  for (size_t i=0; i< 5*4/2; ++i) {
+    CPPUNIT_ASSERT( (args[i].distance >= 0) && (args[i].distance <= 4));
+    if ((args[i].indices.first+1) == args[i].indices.second)
+      CPPUNIT_ASSERT( args[i].bonded );
+  }
+}

src/FunctionApproximation/unittests/ExtractorsUnitTest.hpp

@@ -23 +23 @@
     CPPUNIT_TEST_SUITE( ExtractorsTest) ;
     CPPUNIT_TEST ( gatherAllSymmetricDistancesTest );
+    CPPUNIT_TEST ( gatherAllSymmetricDistances_edgesTest );
     CPPUNIT_TEST_SUITE_END();
 
@@ -29 +30 @@
       void tearDown();
       void gatherAllSymmetricDistancesTest();
+      void gatherAllSymmetricDistances_edgesTest();
 
 private:

src/FunctionApproximation/unittests/Makefile.am

@@ -20 +20 @@
 FUNCTIONAPPROXIMATIONLIBS = \
         libUnitTest.la \
-        ../libMolecuilder.la \
         ../libMolecuilderFunctionApproximation.la \
+        ../libMolecuilderFragmentation.la \
+        ../libMolecuilderFragmentation_getFromKeysetStub.la \
         $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
         ${CodePatterns_LIBS} \