/*
 * 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 "CodePatterns/MemDebug.hpp"

#include <sstream>
#include <utility>
#include <vector>
#include <boost/assign.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"

using namespace boost::assign;

FunctionModel::arguments_t
Extractors::gatherAllDistanceArguments(
    const Fragment::positions_t& positions,
    const Fragment::charges_t& charges,
    const size_t globalid)
{
  FunctionModel::arguments_t result;

  // 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 = positions.begin();//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(
          charges[ std::distance(positions.begin(), firstpositer) ],
          charges[ std::distance(positions.begin(), secondpositer) ]
          );
      arg.indices = std::make_pair(
          std::distance(
              positions.begin(), firstpositer),
          std::distance(
              positions.begin(), secondpositer)
          );
      arg.globalid = globalid;
      result.push_back(arg);
    }
  }

  return result;
}

FunctionModel::arguments_t
Extractors::gatherAllSymmetricDistanceArguments(
    const Fragment::positions_t& positions,
    const Fragment::charges_t& charges,
    const size_t globalid)
{
  FunctionModel::arguments_t result;

  // 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(
          charges[ std::distance(positions.begin(), firstpositer) ],
          charges[ std::distance(positions.begin(), secondpositer) ]
          );
      arg.indices = std::make_pair(
          std::distance(
              positions.begin(), firstpositer),
          std::distance(
              positions.begin(), secondpositer)
          );
      arg.globalid = globalid;
      LOG(3, "DEBUG: Created argument " << arg << ".");
      result.push_back(arg);
    }
  }

  return result;
}

Fragment::positions_t Extractors::_detail::gatherPositionsFromTargets(
    const Fragment::positions_t& positions,
    const Fragment::charges_t& charges,
    const chargeiters_t &targets
    )
{
  Fragment::positions_t filtered_positions;
  for (chargeiters_t::const_iterator firstpairiter = targets.begin();
      firstpairiter != targets.end(); ++firstpairiter) {
    Fragment::positions_t::const_iterator positer = positions.begin();
    const size_t steps = std::distance(charges.begin(), *firstpairiter);
    std::advance(positer, steps);
    filtered_positions.push_back(*positer);
  }
  return filtered_positions;
}

FunctionModel::arguments_t Extractors::_detail::gatherDistancesFromTargets(
    const Fragment::positions_t& positions,
    const Fragment::charges_t& charges,
    const chargeiters_t &targets,
    const size_t globalid
    )
{
  Fragment::positions_t filtered_positions;
  Fragment::charges_t filtered_charges;
  for (chargeiters_t::const_iterator firstpairiter = targets.begin();
      firstpairiter != targets.end(); ++firstpairiter) {
    Fragment::positions_t::const_iterator positer = positions.begin();
    const size_t steps = std::distance(charges.begin(), *firstpairiter);
    std::advance(positer, steps);
    filtered_positions.push_back(*positer);
    filtered_charges.push_back(**firstpairiter);
  }
  return Extractors::gatherAllSymmetricDistanceArguments(
      filtered_positions,
      filtered_charges,
      globalid);
}

Extractors::elementcounts_t 
Extractors::_detail::getElementCounts(
    const Fragment::charges_t elements
  )
{
  elementcounts_t elementcounts;
  for (Fragment::charges_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;
}

Extractors::elementtargets_t 
Extractors::_detail::convertElementcountsToTargets(
    const Fragment::charges_t &charges,
    const elementcounts_t &elementcounts
    )
{
  elementtargets_t elementtargets;
  for (elementcounts_t::const_iterator countiter = elementcounts.begin();
      countiter != elementcounts.end();
      ++countiter) {
    chargeiter_t chargeiter = charges.begin();
    const element_t &element = countiter->first;
    const count_t &count = countiter->second;
    for (count_t i = 0; i < count; ++i) {
      chargeiter_t tempiter = std::find(chargeiter, charges.end(), element);
      if (tempiter != charges.end()) {
        // try to insert new list
        std::pair< elementtargets_t::iterator, bool> inserter =
          elementtargets.insert( std::make_pair( countiter->first, chargeiters_t(1, tempiter)) );
        // if already present, append to it
        if (!inserter.second) {
          inserter.first->second.push_back(tempiter);
        } else { // if created, increase vector's reserve to known size
          inserter.first->second.reserve(countiter->second);
        }
        // search from this element onwards then
        chargeiter = ++tempiter;
      } else {
        ELOG(1, "Could not find desired number " << count << " of element "
            << element << " in fragment with " << charges << ".");
        return Extractors::elementtargets_t();
      }
    }
  }
  return elementtargets;
}

Extractors::elementtargets_t
Extractors::_detail::convertChargesToTargetMap(
    const Fragment::charges_t& charges,
    Fragment::charges_t elements
    )
{
  // place each charge into a map
  elementtargets_t completeelementtargets;
  for (chargeiter_t chargeiter = charges.begin();
      chargeiter != charges.end();
      ++chargeiter) {
    std::pair< elementtargets_t::iterator, bool> inserter =
        completeelementtargets.insert( std::make_pair( *chargeiter, chargeiters_t(1, chargeiter)) );
    // if already present, append to it
    if (!inserter.second) {
      inserter.first->second.push_back(chargeiter);
    }
  }
  // pick out desired charges only
  std::sort(elements.begin(), elements.end());
  Fragment::charges_t::iterator eraseiter =
      std::unique(elements.begin(), elements.end());
  elements.erase(eraseiter, elements.end());
  elementtargets_t elementtargets;
  for(Fragment::charges_t::const_iterator iter = elements.begin();
      iter != elements.end();
      ++iter) {
    elementtargets_t::const_iterator finditer = completeelementtargets.find(*iter);
    ASSERT( finditer != completeelementtargets.end(),
        "Extractors::_detail::convertChargesToTargetMap() - no element "+toString(*iter)+" present?");
    std::pair< elementtargets_t::iterator, bool> inserter =
        elementtargets.insert( std::make_pair( finditer->first, finditer->second) );
    ASSERT( inserter.second,
        "Extractors::_detail::convertChargesToTargetMap() - key twice?");
  }
  return elementtargets;
}


Extractors::chargeiters_t 
Extractors::_detail::realignElementtargets(
    const elementtargets_t &elementtargets,
    const Fragment::charges_t elements,
    const elementcounts_t &elementcounts
    )
{
  chargeiters_t targets;
  elementcounts_t counts; // how many chargeiters of this element have been used
  if (!elements.empty()) {  // skip if no elements given
    targets.reserve(elements.size());
    for (Fragment::charges_t::const_iterator elementiter = elements.begin();
        elementiter != elements.end(); ++elementiter) {
      const element_t &element = *elementiter;
      count_t &count = counts[element]; // if not present, std::map creates instances with default of 0
#ifndef NDEBUG
      {
        elementcounts_t::const_iterator testiter = elementcounts.find(element);
        ASSERT( (testiter != elementcounts.end()) && (count < testiter->second),
            "Extractors::_detail::realignElementTargets() - we want to use more chargeiters for element "
            +toString(element)+" than we counted initially.");
      }
#endif
      elementtargets_t::const_iterator targetiter = elementtargets.find(element);
      if (targetiter != elementtargets.end()) {
        const chargeiters_t &chargeiters = targetiter->second;
        const chargeiter_t &chargeiter = chargeiters[count++];
        targets.push_back(chargeiter);
      }
    }
  }
  return targets;
}

FunctionModel::arguments_t
Extractors::gatherAllDistancesFromFragment(
    const Fragment::positions_t& positions,
    const Fragment::charges_t& charges,
    const Fragment::charges_t elements,
    const size_t globalid
    )
{
  /// The main problem here is that we have to know how many same
  /// elements (but different atoms!) we are required to find. Hence,
  /// we first have to count same elements, then get different targets
  /// for each and then associated them in correct order back again.

  // 0. if no elements given, we return empty arguments
  if (elements.empty())
    return FunctionModel::arguments_t();

  // 1. we have to place each charge into a map as unique chargeiter, i.e. map
  elementtargets_t elementtargets =
      Extractors::_detail::convertChargesToTargetMap(
          charges,
          elements);

  // 2. now we have to combine elementcounts out of elementtargets per desired element
  //    in a combinatorial fashion
  targets_per_combination_t combinations =
      Extractors::_detail::CombineChargesAndTargets(
          elements,
          elementtargets);

  // 3. finally, convert chargeiters into argument list
  FunctionModel::arguments_t args =
      Extractors::_detail::convertTargetsToArguments(
          positions,
          charges,
          combinations,
          globalid);

  return args;
}

Extractors::targets_per_combination_t
Extractors::_detail::CombineChargesAndTargets(
    const Fragment::charges_t& elements,
    const elementtargets_t& elementtargets
    )
{
  // recursively create all correct combinations of targets
  targets_per_combination_t combinations;
  chargeiters_t currenttargets;
  boost::function<void (const chargeiters_t &currenttargets)> addFunction =
      boost::bind(&targets_per_combination_t::push_back,
          boost::ref(combinations),
          _1);
  pickLastElementAsTarget(elements, elementtargets, currenttargets, addFunction);

  return combinations;
}

const Fragment::position_t& getPositionToChargeIter(
    const Fragment::positions_t& positions,
    const Fragment::charges_t& charges,
    const Extractors::chargeiter_t &iter
    )
{
  Fragment::positions_t::const_iterator positer = positions.begin();
  std::advance(positer, std::distance(charges.begin(), iter));
  const Fragment::position_t &position = *positer;
  return position;
}


FunctionModel::arguments_t
Extractors::_detail::convertTargetsToArguments(
    const Fragment::positions_t& positions,
    const Fragment::charges_t& charges,
    const targets_per_combination_t combinations,
    const size_t globalid
    )
{
  FunctionModel::arguments_t args;
  // create arguments from each combination. We cannot use
  // gatherallSymmetricDistanceArguments() because it would not create the
  // correct indices.
  for (targets_per_combination_t::const_iterator iter = combinations.begin();
      iter != combinations.end();
      ++iter) {
    for(chargeiters_t::const_iterator firstiter = iter->begin();
        firstiter != iter->end();
        ++firstiter) {
      const Fragment::position_t &firstpos =
          getPositionToChargeIter(positions, charges, *firstiter);
      const Vector firsttemp(firstpos[0],firstpos[1],firstpos[2]);
      for(chargeiters_t::const_iterator seconditer = firstiter;
          seconditer != iter->end();
          ++seconditer) {
        if (firstiter == seconditer)
          continue;
        const Fragment::position_t &secondpos =
            getPositionToChargeIter(positions, charges, *seconditer);
        const Vector secondtemp(secondpos[0],secondpos[1],secondpos[2]);
        argument_t arg;
        arg.distance = firsttemp.distance(secondtemp);
        arg.indices.first = std::distance(charges.begin(), *firstiter);
        arg.indices.second = std::distance(charges.begin(), *seconditer);
        arg.types.first = **firstiter;
        arg.types.second = **seconditer;
        args.push_back( arg );
      }
    }
  }

  return args;
}

void
Extractors::_detail::pickLastElementAsTarget(
    Fragment::charges_t elements,
    elementtargets_t elementtargets,
    chargeiters_t &currenttargets,
    boost::function<void (const chargeiters_t &currenttargets)> &addFunction
    )
{
  // get last element from charges
  ASSERT( !elements.empty(),
      "Extractors::_detail::pickLastElementAsTarget() - no elements given to pick targets for.");
  const Fragment::charge_t charge = elements.back();
  elements.pop_back();
  elementtargets_t::iterator iter = elementtargets.find(charge);
  if (iter == elementtargets.end())
    return;
  bool NotEmpty = !iter->second.empty();
  while (NotEmpty) {
    // get last target from the vector of chargeiters
    chargeiter_t target = iter->second.back();
    iter->second.pop_back();
    // remove this key if empty
    if (iter->second.empty()) {
      elementtargets.erase(iter);
      NotEmpty = false;
    }
    currenttargets.push_back(target);
    if (elements.empty()) {
      // call add function
      {
        std::stringstream targetstream;
        BOOST_FOREACH( chargeiter_t target, currenttargets ) {
          targetstream << " " << *target;
        }
        LOG(3, "DEBUG: Adding set" << targetstream.str() << ".");
      }
      addFunction(currenttargets);
    } else {
      // if not, call us recursively
      pickLastElementAsTarget(elements, elementtargets, currenttargets, addFunction);
    }
    // pop last in currenset again
    currenttargets.pop_back();
  }
}

Extractors::chargeiters_t
Extractors::_detail::gatherTargetsFromFragment(
    const Fragment::charges_t& charges,
    const Fragment::charges_t elements
    )
{
  /// The main problem here is that we have to know how many same
  /// elements (but different atoms!) we are required to find. Hence,
  /// we first have to count same elements, then get different targets
  /// for each and then associated them in correct order back again.

  // 1. we have to make elements unique with counts, hence convert to map
  elementcounts_t elementcounts =
      Extractors::_detail::getElementCounts(elements);

  // 2. then for each element we need as many targets (chargeiters) as counts
  elementtargets_t elementtargets =
      Extractors::_detail::convertElementcountsToTargets(charges, elementcounts);

  // 3. we go again through elements and use one found target for each count
  // in that order
  chargeiters_t targets =
      Extractors::_detail::realignElementtargets(elementtargets, elements, elementcounts);

#ifndef NDEBUG
  // check all for debugging
  for (chargeiters_t::const_iterator chargeiter = targets.begin();
      chargeiter != targets.end();
      ++chargeiter)
    ASSERT( *chargeiter != charges.end(),
        "Extractors::gatherTargetsFromFragment() - we have not found enough targets?!");
#endif

  return targets;
}

Fragment::positions_t
Extractors::gatherPositionsFromFragment(
    const Fragment::positions_t positions,
    const Fragment::charges_t charges,
    const Fragment::charges_t& elements
    )
{
  // 1.-3. gather correct charge positions
  chargeiters_t targets =
      Extractors::_detail::gatherTargetsFromFragment(charges, elements);
  // 4. convert position_t to Vector
  return Extractors::_detail::gatherPositionsFromTargets(
          positions,
          charges,
          targets);
}

FunctionModel::arguments_t
Extractors::gatherDistancesFromFragment(
    const Fragment::positions_t positions,
    const Fragment::charges_t charges,
    const Fragment::charges_t& elements,
    const size_t globalid
    )
{
  // 1.-3. gather correct charge positions
  chargeiters_t targets =
      Extractors::_detail::gatherTargetsFromFragment(charges, elements);
  // 4. convert position_t to Vector
  return Extractors::_detail::gatherDistancesFromTargets(
          positions,
          charges,
          targets,
          globalid);
}

FunctionModel::list_of_arguments_t Extractors::reorderArgumentsByIncreasingDistance(
    const FunctionModel::list_of_arguments_t &listargs
    )
{
  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;
    FunctionModel::arguments_t sortedargs(args);
    std::sort(sortedargs.begin(), sortedargs.end(), argument_t::DistanceComparator);
    returnargs.push_back(sortedargs);
  }
  return returnargs;
}


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

FunctionModel::list_of_arguments_t Extractors::reorderArgumentsByParticleTypes(
    const FunctionModel::list_of_arguments_t &listargs,
    const ParticleTypes_t &_types
    )
{
  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::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;
          }
          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 << ".");
                break;
              }
            }
          }
        }
        // move along in indices and check bounds
        ++secondindex;
        if (secondindex == indices.end()) {
          ++firstindex;
          if (firstindex != indices.end()-1)
            secondindex = firstindex+1;
        }
      }
      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);
  }

  return returnargs;
}

FunctionModel::list_of_arguments_t Extractors::filterArgumentsByParticleTypes(
    const FunctionModel::arguments_t &args,
    const ParticleTypes_t &_types
    )
{
  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 sortedargs =
      reorderArgumentsByParticleTypes(singlelist_allargs, _types);
  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;
}