SortingPairBuilder.h

Go to the documentation of this file.

#ifndef SORTING_PAIR_BUILDER_H
#define SORTING_PAIR_BUILDER_H

#include <algorithm>

//predicate used to sort index array with std::sort
template<class A>
struct ProxyPred_t
{
    ProxyPred_t(const A &a, unsigned d) : array(a), dim(d) { }

    template<class T>
    bool operator()(const T &a, const T &b)
    { return array[a][dim] < array[b][dim]; }

    unsigned dim;
    const A &array;
};

//simplifies usage by allowing the omission of the template parameter
template<class A>
ProxyPred_t<A> ProxyPred(const A &a, unsigned d)
{
    return ProxyPred_t<A>(a, d);
}

template<class PBase>
class SortingPairBuilder
{
public:
    enum { Dim = PBase::Dim };
    typedef typename PBase::Position_t      Position_t;

    SortingPairBuilder(PBase &p) : particles(p) { }

    template<class Pred, class OP>
    void for_each(const Pred& pred, const OP &op)
    {
        std::size_t size = particles.getLocalNum()+particles.getGhostNum();

        Position_t mean[Dim];
        Position_t variance[Dim];

        //calculate mean position
        std::fill(mean, mean+Dim, 0);
        for(std::size_t i = 0;i<size;++i)
            for(int d = 0;d<Dim;++d)
                mean[d] += particles.R[i][d];

        for(int d = 0;d<Dim;++d)
            mean[d] /= size;

        //calculate variance for each dimension
        std::fill(variance, variance+Dim, 0);
        for(std::size_t i = 0;i<size;++i)
            for(int d = 0;d<Dim;++d)
                variance[d] += (mean[d]-particles.R[i][d])*(mean[d]-particles.R[i][d]);

        int dimension = 0;
        int var = variance[0];
        for(int d = 1;d<Dim;++d)
            if(variance[d]>var)
            {
                dimension = d;
                var = variance[d];
            }

        //sort index array
        std::size_t *indices = new std::size_t[size];
        for(std::size_t i = 0;i<size;++i)
            indices[i] = i;

        std::sort(indices, indices+size, ProxyPred(particles.R, dimension));

        for(std::size_t i = 0;i<size;++i)
            for(std::size_t j = i+1;j<size;++j)
                if(pred(particles.R[indices[i]], particles.R[indices[j]]))
                    op(indices[i], indices[j], particles);
                else if(particles.R[indices[j]][dimension] - particles.R[indices[i]][dimension]
                        > pred.getRange(dimension))
                    break;
    }
private:
    PBase &particles;
};

#endif