/*
* Project: MoleCuilder
* Description: creates and alters molecular systems
* Copyright (C) 2012 University of Bonn. 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 "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.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 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.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 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.indices = std::make_pair(
std::distance(
positions.begin(), firstpositer),
std::distance(
positions.begin(), secondpositer)
);
arg.globalid = globalid;
result.push_back(arg);
}
}
return result;
}
Fragment::positions_t Extractors::_detail::gatherPositionsFromCharges(
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;
}
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 of elements " << count << " in fragment.");
return Extractors::elementtargets_t();
}
}
}
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
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);
ASSERT (targetiter != elementtargets.end(),
"Extractors::_detail::realignElementTargets() - not enough chargeiters for element "
+toString(element)+".");
const chargeiters_t &chargeiters = targetiter->second;
const chargeiter_t &chargeiter = chargeiters[count++];
targets.push_back(chargeiter);
}
return targets;
}
Fragment::positions_t
Extractors::gatherPositionOfTuples(
const Fragment& fragment,
const Fragment::charges_t elements
) {
const Fragment::charges_t charges = fragment.getCharges();
/// 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::gatherDistanceOfTuples() - we have not found enough targets?!");
#endif
// 4. convert position_t to Vector
return Extractors::_detail::gatherPositionsFromCharges(
fragment.getPositions(),
charges,
targets);
}
FunctionModel::arguments_t Extractors::reorderArgumentsByIncreasingDistance(
const FunctionModel::arguments_t &args
)
{
FunctionModel::arguments_t returnargs(args);
std::sort(returnargs.begin(), returnargs.end(), argument_t::DistanceComparator);
return returnargs;
}