RDFCodeCalculator.hpp

Go to the documentation of this file.

/* 
 * RDFCodeCalculator.hpp 
 *
 * This file is part of the Chemical Data Processing Toolkit
 *
 * Copyright (C) 2003 Thomas Seidel <thomas.seidel@univie.ac.at>
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library 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
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; see the file COPYING. If not, write to
 * the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */
 
#ifndef CDPL_DESCR_RDFCODECALCULATOR_HPP
#define CDPL_DESCR_RDFCODECALCULATOR_HPP
 
#include <cstddef>
#include <iterator>
#include <cmath>
#include <functional>
 
#include "CDPL/Math/Matrix.hpp"
#include "CDPL/Math/Vector.hpp"
 
 
namespace CDPL
{
 
    namespace Descr
    {
 
        template <typename T>
        class RDFCodeCalculator
        {
 
          public:
            typedef T EntityType;
 
            typedef std::function<double(const EntityType&, const EntityType&)> EntityPairWeightFunction;
 
            typedef std::function<const Math::Vector3D&(const EntityType&)> Entity3DCoordinatesFunction;
 
            RDFCodeCalculator();
 
            void setSmoothingFactor(double factor);
 
            double getSmoothingFactor() const;
 
            void setScalingFactor(double factor);
 
            double getScalingFactor() const;
 
            void setStartRadius(double start_radius);
 
            double getStartRadius() const;
 
            void setRadiusIncrement(double radius_inc);
 
            double getRadiusIncrement() const;
 
            void setNumSteps(std::size_t num_steps);
 
            std::size_t getNumSteps() const;
 
            void setEntityPairWeightFunction(const EntityPairWeightFunction& func);
 
            void setEntity3DCoordinatesFunction(const Entity3DCoordinatesFunction& func);
 
            void enableDistanceToIntervalCenterRounding(bool enable);
 
            bool distanceToIntervalsCenterRoundingEnabled() const;
 
            template <typename Iter, typename Vec>
            void calculate(Iter beg, Iter end, Vec& rdf_code);
 
          private:
            template <typename Iter>
            void init(Iter beg, Iter end);
 
            double                      smoothingFactor;
            double                      scalingFactor;
            double                      startRadius;
            double                      radiusIncrement;
            std::size_t                 numSteps;
            std::size_t                 numEntities;
            EntityPairWeightFunction    weightFunc;
            Entity3DCoordinatesFunction coordsFunc;
            bool                        distToIntervalCenterRounding;
            Math::DMatrix               weightMatrix;
            Math::DMatrix               distMatrix;
        };
    } // namespace Descr
} // namespace CDPL
 
 
// Implementation
 
template <typename T>
CDPL::Descr::RDFCodeCalculator<T>::RDFCodeCalculator():
    smoothingFactor(1.0), scalingFactor(1.0), startRadius(0.0), radiusIncrement(0.1), numSteps(99),
    distToIntervalCenterRounding(false)
{}
 
template <typename T>
void CDPL::Descr::RDFCodeCalculator<T>::setSmoothingFactor(double factor)
{
    smoothingFactor = factor;
}
 
template <typename T>
void CDPL::Descr::RDFCodeCalculator<T>::setScalingFactor(double factor)
{
    scalingFactor = factor;
}
 
template <typename T>
void CDPL::Descr::RDFCodeCalculator<T>::setStartRadius(double start_radius)
{
    startRadius = start_radius;
}
 
template <typename T>
void CDPL::Descr::RDFCodeCalculator<T>::setRadiusIncrement(double radius_inc)
{
    radiusIncrement = radius_inc;
}
 
template <typename T>
void CDPL::Descr::RDFCodeCalculator<T>::setNumSteps(std::size_t num_steps)
{
    numSteps = num_steps;
}
 
template <typename T>
double CDPL::Descr::RDFCodeCalculator<T>::getSmoothingFactor() const
{
    return smoothingFactor;
}
 
template <typename T>
double CDPL::Descr::RDFCodeCalculator<T>::getScalingFactor() const
{
    return scalingFactor;
}
 
template <typename T>
double CDPL::Descr::RDFCodeCalculator<T>::getStartRadius() const
{
    return startRadius;
}
 
template <typename T>
double CDPL::Descr::RDFCodeCalculator<T>::getRadiusIncrement() const
{
    return radiusIncrement;
}
 
template <typename T>
std::size_t CDPL::Descr::RDFCodeCalculator<T>::getNumSteps() const
{
    return numSteps;
}
 
template <typename T>
void CDPL::Descr::RDFCodeCalculator<T>::setEntityPairWeightFunction(const EntityPairWeightFunction& func)
{
    weightFunc = func;
}
 
template <typename T>
void CDPL::Descr::RDFCodeCalculator<T>::setEntity3DCoordinatesFunction(const Entity3DCoordinatesFunction& func)
{
    coordsFunc = func;
}
 
template <typename T>
void CDPL::Descr::RDFCodeCalculator<T>::enableDistanceToIntervalCenterRounding(bool enable)
{
    distToIntervalCenterRounding = enable;
}
 
template <typename T>
bool CDPL::Descr::RDFCodeCalculator<T>::distanceToIntervalsCenterRoundingEnabled() const
{
    return distToIntervalCenterRounding;
}
 
template <typename T>
template <typename Iter, typename Vec>
void CDPL::Descr::RDFCodeCalculator<T>::calculate(Iter beg, Iter end, Vec& rdf_code)
{
    init(beg, end);
 
    double r = startRadius;
 
    for (std::size_t i = 0; i <= numSteps; i++, r += radiusIncrement) {
        double sum = 0.0;
 
        for (std::size_t j = 0; j < numEntities; j++) {
            for (std::size_t k = j + 1; k < numEntities; k++) {
                double t = r - distMatrix(j, k);
 
                sum += weightMatrix(j, k) * std::exp(-smoothingFactor * t * t);
            }
        }
 
        rdf_code[i] = scalingFactor * sum;
    }
}
 
template <typename T>
template <typename Iter>
void CDPL::Descr::RDFCodeCalculator<T>::init(Iter beg, Iter end)
{
    numEntities = std::distance(beg, end);
 
    weightMatrix.resize(numEntities, numEntities, false);
    distMatrix.resize(numEntities, numEntities, false);
 
    Math::Vector3D entity1_pos;
 
    for (std::size_t i = 0; beg != end; i++) {
        const EntityType& entity1 = *beg;
 
        if (coordsFunc)
            entity1_pos = coordsFunc(entity1);
 
        std::size_t j = i + 1;
 
        for (Iter it = ++beg; it != end; ++it, j++) {
            const EntityType& entity2 = *it;
 
            if (!weightFunc)
                weightMatrix(i, j) = 1.0;
            else
                weightMatrix(i, j) = weightFunc(entity1, entity2);
 
            if (coordsFunc) {
                double dist = length(entity1_pos - coordsFunc(entity2));
 
                if (distToIntervalCenterRounding)
                    dist = std::round((dist - startRadius) / radiusIncrement) * radiusIncrement + startRadius;
 
                distMatrix(i, j) = dist;
 
            } else
                distMatrix(i, j) = 0.0;
        }
    }
}
 
#endif // CDPL_DESCR_RDFCODECALCULATOR_HPP