Changeset 1155ba for src/FunctionApproximation

Timestamp:

Nov 11, 2016, 2:25:36 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:

f45a91

Parents:

f5dbea

git-author:

Frederik Heber <heber@…> (10/08/16 07:56:22)

git-committer:

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

Message:

Converted filterArguments... into faster filter, renamed other to ..ByBindingModel.

• there are now two filters: One filters by list of particles types and is meant for simple pair potentials. The other one filters (and reorders) by the binding model. This takes more computation time but is required for more complex potentials.
• TESTFIX: This allows morse and harmonic pair potentials tests to pass again, removed XFAIL from the corresponding regression tests.

Location:

src/FunctionApproximation

Files:

• Extractors.cpp (modified) (6 diffs)
• Extractors.hpp (modified) (3 diffs)

src/FunctionApproximation/Extractors.cpp

@@ -531 +531 @@
 }
 
-FunctionModel::list_of_arguments_t Extractors::reorderArgumentsByParticleTypes(
-    const FunctionModel::list_of_arguments_t &listargs,
+FunctionModel::list_of_arguments_t Extractors::filterArgumentsByBindingModel(
+    const FunctionModel::arguments_t &args,
     const HomologyGraph &_graph,
     const ParticleTypes_t &_types,
@@ -539 +539 @@
 {
   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 one or no particle types, needs no distances.");
-      continue;
-    }
-    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);
-      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 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);
+
+  // 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;
           }
         }
-        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());
-      }
+        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());
     }
   }
@@ -769 +765 @@
     )
 {
+  // 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
@@ -779 +775 @@
 //  typedef ParticleTypes_t firsttype;
 //  typedef ParticleTypes_t secondtype;
-//  typedef std::map< firsttype, std::map< secondtype, boost::ref(args) > > ArgsLookup_t;
+//  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::arguments_t allargs;
-  allargs.reserve(args.size());
+  FunctionModel::list_of_arguments_t returnargs;
   for (ParticleTypes_t::const_iterator firstiter = _types.begin();
       firstiter != _types.end();
@@ -804 +803 @@
         if ((iter->types.first == *firstiter)
               && (iter->types.second == *seconditer)) {
-          allargs.push_back( *iter );
+          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)) {
-          allargs.push_back( *iter );
+          returnargs.push_back( FunctionModel::arguments_t(1, *iter) );
           iter = availableList.erase(iter);
           LOG(4, "DEBUG: Accepted (flipped) argument.");
@@ -819 +818 @@
     }
   }
-  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.");

src/FunctionApproximation/Extractors.hpp

@@ -103 +103 @@
   }
 
-  /** Reorder the arguments to bring adjacent ones together.
+  /** Filter the arguments to select only these required by the model.
    *
-   * After filtering via particle types arguments related via same indices
-   * must not necessarily be contained in the same bunch. This reordering
-   * is done here, preserving the alignment given in
-   * \sa filterArgumentsByParticleTypes()
+   * \warning this is meant as a faster way of getting the arguments for simple
+   * pair potentials. In any other case, one should use filterArgumentsByBindingModel()
    *
    * \param listargs list of arguments to reorder each
@@ -115 +113 @@
    * \return reordered args
    */
-  FunctionModel::list_of_arguments_t reorderArgumentsByParticleTypes(
-      const FunctionModel::list_of_arguments_t &eachargs,
+  FunctionModel::list_of_arguments_t filterArgumentsByParticleTypes(
+      const FunctionModel::arguments_t &args,
       const HomologyGraph &_graph,
       const ParticleTypes_t &_types,
@@ -122 +120 @@
       );
 
-  /** Filter arguments according to types, allowing multiples.
+  /** Filter and reorder the arguments to bring adjacent ones together.
    *
-   * If particle types is (0,1,2) and three arguments, each with a pair of types,
-   * are given, then the alignment will be: (0,1), (0,2), and (1,2).
+   * We need to find all matching subgraphs (given by \a _bindingmodel) in the
+   * given homology graph (given by \a _graph) of the fragment molecule.
+   * This means filtering down to the desired particle types and then find
+   * all possible matching subgraphs in each of argument lists, \a eachargs.
    *
-   * \param args arguments to reorder
+   * \param listargs list of arguments to filter and order appropriately
    * \param _graph contains binding model of graph
    * \param _types particle type vector
-   * \return filtered list of args
+   * \return reordered args
    */
-  FunctionModel::list_of_arguments_t filterArgumentsByParticleTypes(
+  FunctionModel::list_of_arguments_t filterArgumentsByBindingModel(
       const FunctionModel::arguments_t &args,
       const HomologyGraph &_graph,