/*
* 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 .
*/
/*
* Extractors.cpp
*
* Created on: 15.10.2012
* Author: heber
*/
// include config.h
#ifdef HAVE_CONFIG_H
#include
#endif
#include "CodePatterns/MemDebug.hpp"
#include
#include
#include
#include
#include
#include
#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 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 ¤ttargets,
boost::function &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 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 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(
MapValueIterator(range_t.first),
MapValueIterator(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;
}