dec_optimizer.hpp

#ifndef DEC_OPTIMIZER_HPP
#define DEC_OPTIMIZER_HPP

#include "Grid/iterators/grid_key_dx_iterator_sub.hpp"
#include "Grid/iterators/grid_skin_iterator.hpp"

template <unsigned int dim>
class wavefront : public Box<dim,size_t>
{
public:

    static const int start = 0;

    static const int stop = 1;
};

template <unsigned int dim, typename Graph>
class dec_optimizer
{
    grid_sm<dim,void> gh;

private:

    void expand_one_wf(openfpm::vector<wavefront<dim>> & v_w, std::vector<comb<dim>> & w_comb , size_t d)
    {
        for (size_t j = 0 ; j < dim ; j++)
        {
            v_w.template get<wavefront<dim>::stop>(d)[j] = v_w.template get<wavefront<dim>::stop>(d)[j] + w_comb[d].c[j];
            v_w.template get<wavefront<dim>::start>(d)[j] = v_w.template get<wavefront<dim>::start>(d)[j] + w_comb[d].c[j];
        }
    }


    void adjust_others_wf(openfpm::vector<wavefront<dim>> & v_w,  HyperCube<dim> & hyp, std::vector<comb<dim>> & w_comb, size_t d)
    {
        // expand the intersection of the wavefronts

        std::vector<comb<dim>> q_comb = SubHyperCube<dim,dim-1>::getCombinations_R(w_comb[d],dim-2);

        // Eliminate the w_comb[d] direction

        for (size_t k = 0 ; k < q_comb.size() ; k++)
        {
            for (size_t j = 0 ; j < dim ; j++)
            {
                if (w_comb[d].c[j] != 0)
                {
                    q_comb[k].c[j] = 0;
                }
            }
        }

        // for all the combinations
        for (size_t j = 0 ; j < q_comb.size() ; j++)
        {
            size_t id = hyp.LinId(q_comb[j]);

            // get the combination of the direction d
            bool is_pos = hyp.isPositive(d);

            // is positive, modify the stop point or the starting point

            for (size_t s = 0 ; s < dim ; s++)
            {
                if (is_pos == true)
                {v_w.template get<wavefront<dim>::stop>(id)[s] = v_w.template get<wavefront<dim>::stop>(id)[s] + w_comb[d].c[s];}
                else
                {v_w.template get<wavefront<dim>::start>(id)[s] = v_w.template get<wavefront<dim>::start>(id)[s] + w_comb[d].c[s];}
            }
        }
    }

    template<unsigned int prp> void write_wavefront(Graph & graph,openfpm::vector<wavefront<dim>> & v_w)
    {
        // fill the wall domain with 0

        fill_domain<prp>(graph,gh.getBox(),0);

        // fill the wavefront

        for (int i = 0 ; i < v_w.size() ; i++)
        {
            Box<dim,size_t> box = wavefront<dim>::getBox(v_w.get(i));

            fill_domain<prp>(graph,box,1);
        }
    }

    template<unsigned int p_sub> void fill_domain(Graph & graph,const Box<dim,size_t> & box, long int ids)
    {
        // Create a subgrid iterator
        grid_key_dx_iterator_sub<dim,no_stencil,do_not_print_warning_on_adjustment<dim>> g_sub(gh,box.getKP1(),box.getKP2());

        // iterate through all grid points

        while (g_sub.isNext())
        {
            // get the actual key
            const grid_key_dx<dim> & gk = g_sub.get();

            // get the vertex and set the sub id

            graph.vertex(gh.LinId(gk)).template get<p_sub>() = ids;

            // next subdomain
            ++g_sub;
        }
    }

    template<unsigned int p_sub, unsigned int p_id> void add_to_queue(openfpm::vector<size_t> & domains, openfpm::vector<wavefront<dim>> & v_w, Graph & graph,  std::vector<comb<dim>> & w_comb, long int pr_id, const size_t(& bc)[dim])
    {
        // create a new queue
        openfpm::vector<size_t> domains_new;

        // push in the new queue, the old domains of the queue that are not assigned element

        for (size_t j = 0 ; j < domains.size() ; j++)
        {
            long int gs = graph.vertex(domains.get(j)).template get<p_sub>();
            if (gs < 0)
            {
                // not assigned push it

                domains_new.add(domains.get(j));
            }
        }

        // Create an Hyper-cube
        HyperCube<dim> hyp;

        for (size_t d = 0 ; d < v_w.size() ; d++)
        {
            expand_one_wf(v_w,w_comb,d);
            adjust_others_wf(v_w,hyp,w_comb,d);
        }

        // for each expanded wavefront create a sub-grid iterator and add the sub-domain

        for (size_t d = 0 ; d < v_w.size() ; d++)
        {
            // Create a sub-grid iterator
            grid_key_dx_iterator_sub_bc<dim,no_stencil,do_not_print_warning_on_adjustment<dim>> g_sub(gh,v_w.template get<wavefront<dim>::start>(d),v_w.template get<wavefront<dim>::stop>(d),bc);

            // iterate through all grid points

            while (g_sub.isNext())
            {
                // get the actual key
                const grid_key_dx<dim> & gk = g_sub.get();

                // get the vertex and if does not have a sub-id and is assigned ...
                long int pid = graph.vertex(gh.LinId(gk)).template get<p_sub>();

                // Get the processor id of the sub-sub-domain
                long int pp_id = graph.vertex(gh.LinId(gk)).template get<p_id>();

                // if the sub-sub-domain is not assigned
                if (pid < 0)
                {
                    // ... and we are not processing the full graph
                    if (pr_id != -1)
                    {
                        // ... and the processor id of the sub-sub-domain match the part we are processing, add to the queue

                        if ( pr_id == pp_id)
                            domains_new.add(gh.LinId(gk));
                    }
                    else
                        domains_new.add(gh.LinId(gk));
                }

                ++g_sub;
            }
        }

        // copy the new queue to the old one (it not copied, C++11 move semantic)
        domains.swap(domains_new);
    }

    template <unsigned int p_sub, unsigned int p_id> void expand_from_point(size_t start_p, Graph & graph, Box<dim,size_t> & box, openfpm::vector<wavefront<dim>> & v_w , std::vector<comb<dim>> & w_comb)
    {
        // We assume that Graph is the rapresentation of a cartesian graph
        // this mean that the direction d is at the child d

        // Get the number of wavefronts
        size_t n_wf = w_comb.size();

        // Create an Hyper-cube
        HyperCube<dim> hyp;

        // direction of expansion

        size_t domain_id = graph.vertex(start_p).template get<p_id>();
        bool can_expand = true;

        // while is possible to expand

        while (can_expand)
        {
            // reset can expand
            can_expand = false;

            // for each direction of expansion expand the wavefront

            for (size_t d = 0 ; d < n_wf ; d++)
            {
                // number of processed sub-domain
                size_t n_proc_sub = 0;

                // flag to indicate if the wavefront can expand
                bool w_can_expand = true;

                // Create an iterator of the expanded wavefront
                grid_key_dx<dim> start = grid_key_dx<dim>(v_w.template get<wavefront<dim>::start>(d)) + w_comb[d];
                grid_key_dx<dim> stop = grid_key_dx<dim>(v_w.template get<wavefront<dim>::stop>(d)) + w_comb[d];
                grid_key_dx_iterator_sub<dim,no_stencil,do_not_print_warning_on_adjustment<dim>> it(gh,start,stop);

                // for each sub-domain in the expanded wavefront
                while (it.isNext())
                {
                    // get the wavefront sub-domain id
                    size_t sub_w_e = gh.LinId(it.get());

                    // we get the processor id of the neighborhood sub-domain on direction d
                    // (expanded wavefront)
                    size_t exp_p = graph.vertex(sub_w_e).template get<p_id>();

                    // Check if already assigned
                    long int ass = graph.vertex(sub_w_e).template get<p_sub>();

                    // we check if it is the same processor id and is not assigned
                    w_can_expand &= ((exp_p == domain_id) & (ass < 0));

                    // next domain
                    ++it;

                    // increase the number of processed sub-domain
                    n_proc_sub++;
                }

                // if we did not processed sub-domain, we cannot expand
                w_can_expand &= (n_proc_sub != 0);

                // if you can expand one wavefront we did not reach the end
                can_expand |= w_can_expand;

                // if we can expand the wavefront expand it
                if (w_can_expand == true)
                {
                    // expand the wavefront
                    for (size_t j = 0 ; j < dim ; j++)
                    {
                        v_w.template get<wavefront<dim>::stop>(d)[j] = v_w.template get<wavefront<dim>::stop>(d)[j] + w_comb[d].c[j];
                        v_w.template get<wavefront<dim>::start>(d)[j] = v_w.template get<wavefront<dim>::start>(d)[j] + w_comb[d].c[j];
                    }

                    // expand the intersection of the wavefronts

                    if (dim >= 2)
                    {
                        std::vector<comb<dim>> q_comb = SubHyperCube<dim,dim-1>::getCombinations_R(w_comb[d],dim-2);

                        // Eliminate the w_comb[d] direction

                        for (size_t k = 0 ; k < q_comb.size() ; k++)
                        {
                            for (size_t j = 0 ; j < dim ; j++)
                            {
                                if (w_comb[d].c[j] != 0)
                                {
                                    q_comb[k].c[j] = 0;
                                }
                            }
                        }

                        // for all the combinations
                        for (size_t j = 0 ; j < q_comb.size() ; j++)
                        {
                            size_t id = hyp.LinId(q_comb[j]);

                            // get the combination of the direction d

                            bool is_pos = hyp.isPositive(d);

                            // is positive, modify the stop point or the starting point

                            for (size_t s = 0 ; s < dim ; s++)
                            {
                                if (is_pos == true)
                                {v_w.template get<wavefront<dim>::stop>(id)[s] = v_w.template get<wavefront<dim>::stop>(id)[s] + w_comb[d].c[s];}
                                else
                                {v_w.template get<wavefront<dim>::start>(id)[s] = v_w.template get<wavefront<dim>::start>(id)[s] + w_comb[d].c[s];}
                            }
                        }
                    }
                }
            }
        }

        // get back the hyper-cube produced

        for (size_t i = 0 ; i < dim ; i++)
        {
            // get the index of the wavefront direction
            size_t p_f = hyp.positiveFace(i);
            size_t n_f = hyp.negativeFace(i);

            // set the box
            box.setHigh(i,v_w.template get<wavefront<dim>::stop>(p_f)[i]);
            box.setLow(i,v_w.template get<wavefront<dim>::start>(n_f)[i]);
        }
    }

    void InitializeWavefront(grid_key_dx<dim> & start_p, openfpm::vector<wavefront<dim>> & v_w)
    {
        // Wavefront to initialize

        for (size_t i = 0 ; i < v_w.size() ; i++)
        {
            for (size_t j = 0 ; j < dim ; j++)
            {
                v_w.template get<wavefront<dim>::start>(i)[j] = start_p.get(j);
                v_w.template get<wavefront<dim>::stop>(i)[j] = start_p.get(j);
            }
        }
    }

    template<unsigned int p_id, unsigned int p_sub> grid_key_dx<dim> search_seed(Graph & graph, long int id)
    {
        // if no processor is selected return the first point
        if (id < -1)
        {
            grid_key_dx<dim> key;
            key.zero();

            return key;
        }

        // Create a grid iterator
        grid_key_dx_iterator<dim> g_sub(gh);

        // iterate through all grid points

        while (g_sub.isNext())
        {
            // get the actual key
            const grid_key_dx<dim> & gk = g_sub.get();

            // if the subdomain has the id we are searching stop
            if ((long int)graph.vertex(gh.LinId(gk)).template get<p_id>() == id && graph.vertex(gh.LinId(gk)).template get<p_sub>() == -1)
            {
                return gk;
            }

            ++g_sub;
        }

        // If not found return an invalid key
        grid_key_dx<dim> key;
        key.invalid();

        return key;
    }


    template <unsigned int p_sub, unsigned int p_id> size_t optimize(grid_key_dx<dim> & start_p, Graph & graph, long int pr_id, openfpm::vector<Box<dim,size_t>> & lb, openfpm::vector< openfpm::vector<size_t> > & box_nn_processor , const Ghost<dim,long int> & ghe ,const size_t (& bc)[dim], bool init_sub_id = true, size_t sub_id = 0)
    {
        // queue
        openfpm::vector<size_t> v_q;

        // box list 2
        openfpm::vector< openfpm::vector<size_t> > box_nn_processor2;

        // Create an hyper-cube
        HyperCube<dim> hyp;

        // Get the wavefront combinations
        std::vector<comb<dim>> w_comb = hyp.getCombinations_R(dim-1);

        // wavefronts
        openfpm::vector<wavefront<dim>> v_w(w_comb.size());

        // fill the sub decomposition with negative number

        if (init_sub_id == true)
            fill_domain<p_sub>(graph,gh.getBox(),-1);

        // push the first domain
        v_q.add(gh.LinId(start_p));

        while (v_q.size() != 0)
        {
            // Box
            Box<dim,size_t> box;

            // Get the grid_key position from the linearized id
            start_p = gh.InvLinId(v_q.get(0));

            // Initialize the wavefronts from the domain start_p
            InitializeWavefront(start_p,v_w);

            // Create the biggest box containing the domain
            expand_from_point<p_sub,p_id>(v_q.get(0),graph,box,v_w,w_comb);

            // Add the created box to the list of boxes
            lb.add(box);

            // fill the domain
            fill_domain<p_sub>(graph,box,sub_id);

            // add the surrounding sub-domain to the queue
            add_to_queue<p_sub,p_id>(v_q,v_w,graph,w_comb,pr_id,bc);

            // increment the sub_id
            sub_id++;
        }

        return sub_id;
    }

    template<unsigned int p_id> void construct_box_nn_processor(Graph & graph, openfpm::vector< openfpm::vector<size_t> > & box_nn_processor, const openfpm::vector<Box<dim,size_t>> & subs, const Ghost<dim,long int> & ghe, const size_t (& bc)[dim], long int pr_id)
    {
        std::unordered_map<size_t,size_t> map;

        for (size_t i = 0 ; i < subs.size() ; i++)
        {
            map.clear();
            Box<dim,size_t> sub = subs.get(i);
            sub.enlarge(ghe);

            grid_skin_iterator_bc<dim> gsi(gh,subs.get(i),sub,bc);

            while (gsi.isNext())
            {
                auto key = gsi.get();

                size_t pp_id = graph.vertex(gh.LinId(key)).template get<p_id>();
                if (pr_id != (long int)pp_id)
                    map[pp_id] = pp_id;

                ++gsi;
            }

            // Add the keys to box_nn_processors

            box_nn_processor.add();
            for ( auto it = map.begin(); it != map.end(); ++it )
            {
                box_nn_processor.last().add(it->first);
            }
        }
    }

public:

    dec_optimizer(Graph & g, const size_t (& sz)[dim])
    :gh(sz)
    {
        // The graph g is suppose to represent a cartesian grid
        // No check is performed on g
    }

    template <unsigned int p_sub, unsigned int p_id> void optimize(grid_key_dx<dim> & start_p, Graph & graph, const Ghost<dim,long int> & ghe , const size_t (& bc)[dim])
    {
        // temporal vector
        openfpm::vector<Box<dim,size_t>> tmp;

        // temporal vector
        openfpm::vector< openfpm::vector<size_t> > box_nn_processor;

        // optimize
        optimize<p_sub,p_id>(start_p,graph,-1,tmp, box_nn_processor,ghe,bc);
    }

    template <unsigned int p_sub, unsigned int p_id> void optimize(Graph & graph, long int pr_id, openfpm::vector<Box<dim,size_t>> & lb, openfpm::vector< openfpm::vector<size_t> > & box_nn_processor, const Ghost<dim,long int> & ghe, const size_t (& bc)[dim])
    {
        grid_key_dx<dim> key_seed;
        key_seed.zero();

        // if processor is -1 call optimize with -1 to do on all processors and exit
        if (pr_id == -1)
        {
            optimize<p_sub,p_id>(key_seed,graph,pr_id,lb,box_nn_processor,ghe,bc);

            // Construct box box_nn_processor from the constructed domain
            construct_box_nn_processor<p_id>(graph,box_nn_processor,lb,ghe,bc,pr_id);

            return;
        }

        size_t sub_id = 0;

        // fill the sub decomposition with negative number
        fill_domain<p_sub>(graph,gh.getBox(),-1);

        key_seed = search_seed<p_id,p_sub>(graph,pr_id);

        while (key_seed.isValid())
        {
            // optimize
            sub_id = optimize<p_sub,p_id>(key_seed,graph,pr_id,lb,box_nn_processor,ghe,bc,false,sub_id);

            // new seed
            key_seed = search_seed<p_id,p_sub>(graph,pr_id);
        }

        // Construct box box_nn_processor from the constructed domain
        construct_box_nn_processor<p_id>(graph,box_nn_processor,lb,ghe,bc,pr_id);
    }
};

#endif