ORB.hpp

/*
 * ORB.hpp
 *
 *  Created on: Mar 13, 2015
 *      Author: Pietro Incardona
 */
 
#ifndef ORB_HPP_
#define ORB_HPP_
 
#include "util/mathutil.hpp"
 
template<unsigned int dim, typename ORB>
struct bisect_unroll
{
    ORB & orb;
 
    bisect_unroll(ORB & orb)
    :orb(orb)
    {};
 
    template<typename T>
    void operator()(T& t)
    {
        orb.template bisect<T::value>();
    }
};
 
template<unsigned int dim, unsigned int i, typename ORB, class Enable=void>
struct do_when_dim_gr_i
{
    static void bisect_loop(bisect_unroll<dim,ORB> & bu)
    {
    }
};
 
template<unsigned int dim, unsigned int i, typename ORB>
struct do_when_dim_gr_i<dim,i,ORB,typename boost::enable_if< boost::mpl::bool_<(i < dim)> >::type>
{
    static void bisect_loop(bisect_unroll<dim,ORB> & bu, typename boost::enable_if< boost::mpl::bool_<(i < dim)> >::type* dummy = 0)
    {
        boost::mpl::for_each< boost::mpl::range_c<int,0,i> >(bu);
    }
};
 
 
template<typename T> class ORB_node : public aggregate<T>
{
public:
 
    static const unsigned int CM = 0;
};
 
template<unsigned int dim, typename T, typename loc_wg=openfpm::vector<float>, typename loc_pos=openfpm::vector<Point<dim,T>> , typename Box=Box<dim,T>, template<typename,typename> class Tree=Graph_CSR_s>
class ORB
{
    // Virtual cluster
    Vcluster<> & v_cl;
 
    // particle coordinate accumulator
    openfpm::vector<T> cm;
    // number of elements summed into cm
    openfpm::vector<size_t> cm_cnt;
 
    // Local particle vector
    loc_pos & lp;
 
    // Label id of the particle, the label identify where the "local" particles
    // are in the local decomposition
    openfpm::vector<size_t> lp_lbl;
 
    size_t dim_div;
 
    // Structure that store the orthogonal bisection tree
    Tree<ORB_node<T>,no_edge> grp;
 
    template<unsigned int dir> void local_cm(size_t start)
    {
        typedef Point<dim,T> s;
 
        // reset the counters and accumulators
        cm.fill(0);
        cm_cnt.fill(0);
 
        // Calculate the local CM
        auto it = lp.getIterator();
 
        while (it.isNext())
        {
            // vector key
            auto key =  it.get();
 
            size_t lbl = lp_lbl.get(key) - start;
 
            // add the particle coordinate to the CM accumulator
            cm.get(lbl) += lp.template get<s::x>(key)[dir];
 
            // increase the counter
            cm_cnt.get(lbl)++;
 
            ++it;
        }
    }
 
    template<unsigned int dir> inline void local_label()
    {
        typedef Point<dim,T> s;
 
        // direction of the previous split
        const size_t dir_cm = (dir == 0)?(dim-1):(dir-1);
 
        // if it is the first iteration we just have to create the labels array with zero
        if (grp.getNVertex() == 1)
        {lp_lbl.resize(lp.size());lp_lbl.fill(0); return;}
 
        // we check where (the node) the particles live
 
        auto p_it = lp.getIterator();
 
        while(p_it.isNext())
        {
            auto key = p_it.get();
 
            // Get the label
            size_t lbl = lp_lbl.get(key);
 
            // each node has two child's (FOR SURE)
            // No leaf are processed here
            size_t n1 = grp.getChild(lbl,0);
            size_t n2 = grp.getChild(lbl,1);
 
            // get the splitting center of mass
            T cm = grp.template vertex_p<ORB_node<T>::CM>(lbl);
 
            // if the particle n-coordinate is smaller than the CM is inside the child n1
            // otherwise is inside the child n2
 
            if (lp.template get<s::x>(key)[dir_cm] < cm)
            {lp_lbl.get(key) = n1;}
            else
            {lp_lbl.get(key) = n2;}
 
            ++p_it;
        }
    }
 
    template<unsigned int dir> size_t bisect()
    {
        //
        size_t start = grp.getNVertex();
 
        // first label the local particles
        local_label<dir>();
 
        // Index from where start the first leaf
        size_t n_node = (start + 1) / 2;
 
        // Calculate the local cm
        local_cm<dir>(start - n_node);
 
        // reduce the local cm and cm_cnt
        v_cl.sum(cm);
        v_cl.sum(cm_cnt);
        v_cl.execute();
 
        // set the CM for the previous leaf (they are not anymore leaf)
 
        for (size_t i = 0 ; i < n_node ; i++)
        {
            grp.template vertex_p<ORB_node<T>::CM>(start-n_node+i) = cm.get(i) / cm_cnt.get(i);
        }
 
        // append on the 2^(n-1) previous end node to the 2^n leaf
 
        for (size_t i = start - n_node, s = 0 ; i < start ; i++, s++)
        {
            // Add 2 leaf nodes and connect them with the node
            grp.addVertex();
            grp.addVertex();
            grp.template addEdge(i,start+2*s);
            grp.template addEdge(i,start+2*s+1);
        }
 
        return 2*n_node;
    }
 
    // define the friend class bisect_unroll
 
    friend class bisect_unroll<dim,ORB<dim,T,loc_wg,loc_pos,Box,Tree>>;
 
public:
 
    ORB(Box dom, size_t n_sub, loc_pos & lp)
    :v_cl(create_vcluster()),lp(lp)
    {
        typedef ORB<dim,T,loc_wg,loc_pos,Box,Tree> ORB_class;
 
        dim_div = 0;
 
        n_sub = openfpm::math::round_big_2(n_sub);
        size_t nsub = log2(n_sub);
 
        // number of center or mass needed
        cm.resize(2 << nsub);
        cm_cnt.resize(2 << nsub);
 
        // Every time we divide from 0 to dim-1, calculate how many cycle from 0 to dim-1 we have
        size_t dim_cycle = nsub / dim;
 
        // Create the root node in the graph
        grp.addVertex();
 
        // unroll bisection cycle
        bisect_unroll<dim,ORB_class> bu(*this);
        for (size_t i = 0 ; i < dim_cycle ; i++)
        {
            boost::mpl::for_each< boost::mpl::range_c<int,0,dim> >(bu);
            // bu is recreated several time internaly
        }
 
        // calculate and execute the remaining cycles
        switch(nsub - dim_cycle * dim)
        {
        case 1:
            do_when_dim_gr_i<dim,1,ORB_class>::bisect_loop(bu);
            break;
        case 2:
            do_when_dim_gr_i<dim,2,ORB_class>::bisect_loop(bu);
            break;
        case 3:
            do_when_dim_gr_i<dim,3,ORB_class>::bisect_loop(bu);
            break;
        case 4:
            do_when_dim_gr_i<dim,4,ORB_class>::bisect_loop(bu);
            break;
        case 5:
            do_when_dim_gr_i<dim,5,ORB_class>::bisect_loop(bu);
            break;
        case 6:
            do_when_dim_gr_i<dim,6,ORB_class>::bisect_loop(bu);
            break;
        case 7:
            do_when_dim_gr_i<dim,7,ORB_class>::bisect_loop(bu);
            break;
        default:
            // To extend on higher dimension create other cases or create runtime
            // version of local_cm and bisect
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " ORB is not working for dimension bigger than 8";
 
        }
    }
};
 
 
#endif /* ORB_HPP_ */