CartDecomposition.hpp

/*
 * CartDecomposition.hpp
 *
 *  Created on: Oct 07, 2015
 *      Author: Pietro Incardona, Antonio Leo
 */
 
#ifndef CARTDECOMPOSITION_HPP
#define CARTDECOMPOSITION_HPP
 
#include "config.h"
#include <cmath>
#include "VCluster/VCluster.hpp"
#include "Graph/CartesianGraphFactory.hpp"
#include "Decomposition.hpp"
#include "Vector/map_vector.hpp"
#include <vector>
#include <initializer_list>
#include "SubdomainGraphNodes.hpp"
#include "dec_optimizer.hpp"
#include "Space/Shape/Box.hpp"
#include "Space/Shape/Point.hpp"
#include "NN/CellList/CellDecomposer.hpp"
#include <unordered_map>
#include "NN/CellList/CellList.hpp"
#include "Space/Ghost.hpp"
#include "common.hpp"
#include "ie_loc_ghost.hpp"
#include "ie_ghost.hpp"
#include "nn_processor.hpp"
#include "GraphMLWriter/GraphMLWriter.hpp"
#include "Distribution/ParMetisDistribution.hpp"
#include "Distribution/DistParMetisDistribution.hpp"
#include "Distribution/MetisDistribution.hpp"
#include "DLB/DLB.hpp"
#include "util/se_util.hpp"
#include "util/mathutil.hpp"
#include "CartDecomposition_ext.hpp"
#include "data_type/aggregate.hpp"
#include "Domain_NN_calculator_cart.hpp"
#include "cuda/CartDecomposition_gpu.cuh"
#include "Domain_icells_cart.hpp"
 
#define CARTDEC_ERROR 2000lu
 
enum dec_options
{
    DEC_NONE = 0,
    DEC_SKIP_ICELL = 1
};
 
template<unsigned int dim> static void nsub_to_div2(size_t (& div)[dim], size_t n_sub, size_t dim_r)
{
    for (size_t i = 0; i < dim; i++)
    {
        if (i < dim_r)
        {div[i] = openfpm::math::round_big_2(pow(n_sub, 1.0 / dim_r));}
        else
        {div[i] = 1;}
    }
}
 
template<unsigned int dim> static void nsub_to_div(size_t (& div)[dim], size_t n_sub, size_t dim_r)
{
    for (size_t i = 0; i < dim; i++)
    {
        if (i < dim_r)
        {div[i] = std::floor(pow(n_sub, 1.0 / dim_r));}
        else
        {div[i] = 1;}
    }
}
 
#define COMPUTE_SKIN_SUB 1
 
template<unsigned int dim, typename T, typename Memory, template <typename> class layout_base, typename Distribution>
class CartDecomposition: public ie_loc_ghost<dim, T,layout_base, Memory>,
                         public nn_prcs<dim, T,layout_base,Memory>,
                         public ie_ghost<dim,T,Memory,layout_base>,
                         public domain_nn_calculator_cart<dim>,
                         public domain_icell_calculator<dim,T,layout_base,Memory>
{
public:
 
    typedef T domain_type;
 
    // typedef Box<dim, T> Box_;
 
    typedef CartDecomposition<dim,T,Memory,layout_base,Distribution> base_type;
 
    typedef CartDecomposition_ext<dim,T,Memory,layout_base,Distribution> extended_type;
 
    typedef ie_loc_ghost<dim, T,layout_base, Memory> ie_loc_ghost_type;
 
    typedef nn_prcs<dim, T,layout_base,Memory> nn_prcs_type;
 
    typedef ie_ghost<dim,T,Memory,layout_base> ie_ghost_type;
 
protected:
 
    bool host_dev_transfer = false;
 
    bool commCostSet = false;
 
    typedef typename openfpm::vector<Box<dim, T>,
            Memory,
            memory_traits_lin,
            openfpm::vector_grow_policy_default,
            openfpm::vect_isel<Box<dim, T>>::value>::access_key acc_key;
 
    openfpm::vector<Box<dim, T>,Memory,layout_base> sub_domains;
 
    mutable openfpm::vector<Box_map<dim, T>,Memory,layout_base> sub_domains_global;
 
    openfpm::vector<openfpm::vector<long unsigned int> > box_nn_processor;
 
    CellList<dim,T,Mem_fast<Memory,int>,shift<dim,T>> fine_s;
 
    grid_sm<dim, void> gr;
 
    grid_sm<dim, void> gr_dist;
 
    CellDecomposer_sm<dim, T, shift<dim,T>> cd;
 
    Box<dim,T> domain;
 
    T spacing[dim];
 
    size_t magn[dim];
 
    Vcluster<> & v_cl;
 
    Distribution dist;
 
    Box<dim,T> bbox;
 
    long int ref_cnt;
 
    Ghost<dim,T> ghost;
 
    size_t bc[dim];
 
    Box<dim,size_t> proc_box;
 
    openfpm::vector<Box<dim, size_t>> loc_boxes;
 
    template<typename Memory_bx> Box<dim,T> convertDecBoxIntoSubDomain(encapc<1,Box<dim,size_t>,Memory_bx> loc_box)
    {
        // A point with all coordinate to one
        size_t one[dim];
        for (size_t i = 0 ; i < dim ; i++)  {one[i] = 1;}
 
        Box<dim, size_t> sub_dc = loc_box;
        Box<dim, size_t> sub_dce = sub_dc;
        sub_dce.expand(one);
        sub_dce.mul(magn);
 
        // shrink by one
        for (size_t i = 0 ; i < dim ; i++)
        {
            loc_box.template get<Box<dim,T>::p1>()[i] = sub_dce.getLow(i);
            loc_box.template get<Box<dim,T>::p2>()[i] = sub_dce.getHigh(i) - 1;
        }
 
        Box<dim, T> sub_d(sub_dce);
        sub_d.mul(spacing);
        sub_d += domain.getP1();
 
        // we add the
 
        // Fixing sub-domains to cover all the domain
 
        // Fixing sub_d
        // if (loc_box) is at the boundary we have to ensure that the box span the full
        // domain (avoiding rounding off error)
        for (size_t i = 0; i < dim; i++)
        {
            if (sub_dc.getHigh(i) == gr.size(i) - 1)
            {sub_d.setHigh(i, domain.getHigh(i));}
 
            if (sub_dc.getLow(i) == 0)
            {sub_d.setLow(i,domain.getLow(i));}
        }
 
        return sub_d;
    }
 
    void collect_all_sub_domains(openfpm::vector<Box_map<dim,T>,Memory,layout_base> & sub_domains_global)
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif
 
        sub_domains_global.clear();
        openfpm::vector<Box_map<dim,T>,Memory,layout_base> bm;
 
        for (size_t i = 0 ; i < sub_domains.size() ; i++)
        {
            bm.add();
 
            bm.template get<0>(bm.size()-1) = Box<dim,T>(sub_domains.get(i));
            bm.template get<1>(bm.size()-1) = v_cl.rank();
        }
 
        v_cl.SGather<decltype(bm),decltype(sub_domains_global),layout_base>(bm,sub_domains_global,0);
 
        size_t size = sub_domains_global.size();
 
        v_cl.max(size);
        v_cl.execute();
 
        sub_domains_global.resize(size);
 
        v_cl.Bcast(sub_domains_global,0);
        v_cl.execute();
    }
 
public:
 
    void initialize_fine_s(const Box<dim,T> & domain)
    {
        fine_s.clear();
        size_t div_g[dim];
 
        // We reduce the size of the cells by a factor 8 in 3d 4 in 2d
        for (size_t i = 0 ; i < dim ; i++)
        {div_g[i] = (gr.size(i) == 1)?1:gr.size(i)/2;}
 
        fine_s.Initialize(domain,div_g);
    }
 
    void construct_fine_s()
    {
        collect_all_sub_domains(sub_domains_global);
 
        // now draw all sub-domains in fine-s
 
        for (size_t i = 0 ; i < sub_domains_global.size() ; i++)
        {
 
            // get the cells this box span
            const grid_key_dx<dim> p1 = fine_s.getCellGrid_me(sub_domains_global.template get<0>(i).getP1());
            const grid_key_dx<dim> p2 = fine_s.getCellGrid_pe(sub_domains_global.template get<0>(i).getP2());
 
            // Get the grid and the sub-iterator
            auto & gi = fine_s.getGrid();
            grid_key_dx_iterator_sub<dim> g_sub(gi,p1,p2);
 
            // add the box-id to the cell list
            while (g_sub.isNext())
            {
                auto key = g_sub.get();
                fine_s.addCell(gi.LinId(key),i);
 
                ++g_sub;
            }
        }
 
        host_dev_transfer = false;
    }
 
    void createSubdomains(Vcluster<> & v_cl, const size_t (& bc)[dim], size_t opt = 0)
    {
        int p_id = v_cl.getProcessUnitID();
 
        // Calculate the total number of box and and the spacing
        // on each direction
        // Get the box containing the domain
        Box<dim, T> bs = domain.getBox();
 
        for (unsigned int i = 0; i < dim; i++)
        {
            // Calculate the spacing
            spacing[i] = (bs.getHigh(i) - bs.getLow(i)) / gr.size(i);
        }
 
        // fill the structure that store the processor id for each sub-domain
        initialize_fine_s(domain);
 
        // Optimize the decomposition creating bigger spaces
        // And reducing Ghost over-stress
        dec_optimizer<dim, Graph_CSR<nm_v<dim>, nm_e>> d_o(dist.getGraph(), gr_dist.getSize());
 
        // Ghost
        Ghost<dim,long int> ghostExtended;
 
        // Set the ghost
        for (size_t i = 0 ; i < dim ; i++)
        {
            ghostExtended.setLow(i,static_cast<long int>(ghost.getLow(i)/spacing[i]) - 1);
            ghostExtended.setHigh(i,static_cast<long int>(ghost.getHigh(i)/spacing[i]) + 1);
        }
 
        // optimize the decomposition or merge sub-sub-domain
        d_o.template optimize<nm_v_sub_id, nm_v_proc_id>(dist.getGraph(), p_id, loc_boxes, box_nn_processor,ghostExtended,bc);
 
        // Initialize
        if (loc_boxes.size() > 0)
        {
            bbox = convertDecBoxIntoSubDomain(loc_boxes.get(0));
            proc_box = loc_boxes.get(0);
            sub_domains.add(bbox);
        }
        else
        {
            // invalidate all the boxes
            for (size_t i = 0 ; i < dim ; i++)
            {
                proc_box.setLow(i,0.0);
                proc_box.setHigh(i,0);
 
                bbox.setLow(i,0.0);
                bbox.setHigh(i,0);
            }
        }
 
        // convert into sub-domain
        for (size_t s = 1; s < loc_boxes.size(); s++)
        {
            Box<dim,T> sub_d = convertDecBoxIntoSubDomain(loc_boxes.get(s));
 
            // add the sub-domain
            sub_domains.add(sub_d);
 
            // Calculate the bound box
            bbox.enclose(sub_d);
            proc_box.enclose(loc_boxes.get(s));
        }
 
        nn_prcs<dim,T,layout_base,Memory>::create(box_nn_processor, sub_domains);
        nn_prcs<dim,T,layout_base,Memory>::applyBC(domain,ghost,bc);
 
        // fill fine_s structure
        // fine_s structure contain the processor id for each sub-sub-domain
        // with sub-sub-domain we mean the sub-domain decomposition before
        // running dec_optimizer (before merging sub-domains)
 
 
        construct_fine_s();
 
        Initialize_geo_cell_lists();
    }
 
    void Initialize_geo_cell_lists()
    {
        // Get the processor bounding Box
        Box<dim,T> bound = getProcessorBounds();
 
        // Check if the box is valid
        if (bound.isValidN() == true)
        {
            // calculate the sub-divisions
            size_t div[dim];
            for (size_t i = 0; i < dim; i++)
            {div[i] = (size_t) ((bound.getHigh(i) - bound.getLow(i)) / cd.getCellBox().getP2()[i]);}
 
            // Initialize the geo_cell structure
            ie_ghost<dim,T,Memory,layout_base>::Initialize_geo_cell(bound,div);
 
            // Initialize shift vectors
            ie_ghost<dim,T,Memory,layout_base>::generateShiftVectors(domain,bc);
        }
    }
 
    void computeCommunicationAndMigrationCosts(size_t ts)
    {
        float migration = 0;
 
        Box<dim, T> cellBox = cd.getCellBox();
        float b_s = static_cast<float>(cellBox.getHigh(0));
        float gh_s = static_cast<float>(ghost.getHigh(0));
 
        // compute the gh_area for 2 dim case
        float gh_v = (gh_s * b_s);
 
        // multiply for sub-sub-domain side for each domain
        for (size_t i = 2; i < dim; i++)
        {
            /* coverity[dead_error_line] */
            gh_v *= b_s;
        }
 
        size_t norm = (size_t) (1.0 / gh_v);
 
        migration = pow(b_s, dim);
 
        size_t prev = 0;
 
        for (size_t i = 0; i < dist.getNSubSubDomains(); i++)
        {
            dist.setMigrationCost(i, norm * migration /* * dist.getSubSubDomainComputationCost(i)*/ );
 
            for (size_t s = 0; s < dist.getNSubSubDomainNeighbors(i); s++)
            {
                // We have to remove dist.getSubSubDomainComputationCost(i) otherwise the graph is
                // not directed
                dist.setCommunicationCost(i, s, 1   *  ts);
            }
            prev += dist.getNSubSubDomainNeighbors(i);
        }
 
        commCostSet = true;
    }
 
    void CreateSubspaces()
    {
        // Create a grid where each point is a space
        grid_sm<dim, void> g(div);
 
        // create a grid_key_dx iterator
        grid_key_dx_iterator<dim> gk_it(g);
 
        // Divide the space into subspaces
        while (gk_it.isNext())
        {
            grid_key_dx<dim> key = gk_it.get();
 
            Box<dim, T> tmp;
 
            for (int i = 0; i < dim; i++)
            {
                tmp.setHigh(i, (key.get(i) + 1) * spacing[i]);
                tmp.setLow(i, key.get(i) * spacing[i]);
            }
 
            sub_domains.add(tmp);
 
            // Next sub-domain
            ++gk_it;
        }
    }
 
 
    void calculateGhostBoxes()
    {
        // Intersect all the local sub-domains with the sub-domains of the contiguous processors
 
        // create the internal structures that store ghost information
        ie_ghost<dim, T,Memory,layout_base>::create_box_nn_processor_ext(v_cl, ghost, sub_domains, box_nn_processor, *this);
        ie_ghost<dim, T,Memory,layout_base>::create_box_nn_processor_int(v_cl, ghost, sub_domains, box_nn_processor, *this);
 
        ie_loc_ghost<dim,T,layout_base,Memory>::create(sub_domains,domain,ghost,bc);
 
        // Ghost box information must be re-offloaded
        host_dev_transfer = false;
        ie_ghost<dim, T,Memory,layout_base>::reset_host_dev_transfer();
    }
 
 
public:
 
    static constexpr int dims = dim;
 
    typedef T stype;
 
    void incRef()
    {ref_cnt++;}
 
    void decRef()
    {ref_cnt--;}
 
    long int ref()
    {
        return ref_cnt;
    }
 
    CartDecomposition(Vcluster<> & v_cl)
    :nn_prcs<dim, T, layout_base,Memory>(v_cl), v_cl(v_cl), dist(v_cl),ref_cnt(0)
    {
        // Reset the box to zero
        bbox.zero();
    }
 
    CartDecomposition(const CartDecomposition<dim,T,Memory,layout_base,Distribution> & cart)
    :nn_prcs<dim,T,layout_base,Memory>(cart.v_cl),v_cl(cart.v_cl),dist(v_cl),ref_cnt(0)
    {
        this->operator=(cart);
    }
 
    CartDecomposition(CartDecomposition<dim,T,Memory,layout_base,Distribution> && cart)
    :nn_prcs<dim,T,layout_base,Memory>(cart.v_cl),v_cl(cart.v_cl),dist(v_cl),ref_cnt(0)
    {
        this->operator=(cart);
    }
 
    ~CartDecomposition()
    {
    }
 
    class box_id
    {
    public:
        inline static size_t id(p_box<dim, T> & p, size_t b_id)
        {
            return b_id;
        }
    };
 
    class processor_id
    {
    public:
        template<typename encap_type> inline static size_t id(const encap_type & p, size_t b_id)
        {
            return p.template get<proc_>();
        }
    };
 
    class lc_processor_id
    {
    public:
        template<typename encap_type> inline static size_t id(const encap_type & p, size_t b_id)
        {
            return p.template get<lc_proc_>();
        }
    };
 
    class shift_id
    {
    public:
        template<typename encap_type> inline static size_t id(const encap_type & p, size_t b_id)
        {
            return p.template get<shift_id_>();
        }
    };
 
    void applyPointBC(float (& pt)[dim]) const
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (bc[i] == PERIODIC)
            {pt[i] = openfpm::math::periodic_l(pt[i],domain.getHigh(i),domain.getLow(i));}
        }
    }
 
    void applyPointBC(Point<dim,T> & pt) const
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (bc[i] == PERIODIC)
            {pt.get(i) = openfpm::math::periodic_l(pt.get(i),domain.getHigh(i),domain.getLow(i));}
        }
    }
 
    template<typename Mem> void applyPointBC(encapc<1,Point<dim,T>,Mem> && pt) const
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (bc[i] == PERIODIC)
            {pt.template get<0>()[i] = openfpm::math::periodic_l(pt.template get<0>()[i],domain.getHigh(i),domain.getLow(i));}
        }
    }
 
    CartDecomposition<dim,T,Memory,layout_base,Distribution> duplicate(const Ghost<dim,T> & g) const
    {
        CartDecomposition<dim,T,Memory,layout_base,Distribution> cart(v_cl);
 
        cart.box_nn_processor = box_nn_processor;
        cart.sub_domains = sub_domains;
        cart.fine_s = fine_s;
 
        cart.gr = gr;
        cart.cd = cd;
        cart.domain = domain;
        for (size_t i = 0 ; i < dim ; i++)
        {
            cart.spacing[i] = spacing[i];
            cart.magn[i] = magn[i];
        };
 
        cart.bbox = bbox;
        cart.ghost = g;
 
        cart.dist = dist;
 
        for (size_t i = 0 ; i < dim ; i++)
            cart.bc[i] = bc[i];
 
        (static_cast<nn_prcs<dim,T,layout_base,Memory> &>(cart)).create(box_nn_processor, sub_domains);
        (static_cast<nn_prcs<dim,T,layout_base,Memory> &>(cart)).applyBC(domain,ghost,bc);
 
        cart.Initialize_geo_cell_lists();
        cart.calculateGhostBoxes();
 
        cart.collect_all_sub_domains(cart.sub_domains_global);
 
        return cart;
    }
 
    CartDecomposition<dim,T,Memory,layout_base,Distribution> duplicate() const
    {
        CartDecomposition<dim,T,Memory,layout_base,Distribution> cart(v_cl);
 
        (static_cast<ie_loc_ghost<dim,T, layout_base,Memory>*>(&cart))->operator=(static_cast<ie_loc_ghost<dim,T,layout_base,Memory>>(*this));
        (static_cast<nn_prcs<dim,T,layout_base,Memory>*>(&cart))->operator=(static_cast<nn_prcs<dim,T,layout_base,Memory>>(*this));
        (static_cast<ie_ghost<dim,T,Memory,layout_base>*>(&cart))->operator=(static_cast<ie_ghost<dim,T,Memory,layout_base>>(*this));
 
        cart.sub_domains = sub_domains;
        cart.box_nn_processor = box_nn_processor;
        cart.fine_s = fine_s;
        cart.gr = gr;
        cart.gr_dist = gr_dist;
        cart.dist = dist;
        cart.commCostSet = commCostSet;
        cart.cd = cd;
        cart.domain = domain;
        cart.sub_domains_global = sub_domains_global;
        for (size_t i = 0 ; i < dim ; i++)
        {
            cart.spacing[i] = spacing[i];
            cart.magn[i] = magn[i];
        };
 
        cart.ghost = ghost;
 
        cart.bbox = bbox;
 
        for (size_t i = 0 ; i < dim ; i++)
        {cart.bc[i] = this->bc[i];}
 
        return cart;
    }
 
    template<typename Memory2, template <typename> class layout_base2>
    CartDecomposition<dim,T,Memory2,layout_base2,Distribution> duplicate_convert() const
    {
        CartDecomposition<dim,T> cart(v_cl);
 
        (static_cast<ie_loc_ghost<dim,T,layout_base2,Memory2>*>(&cart))->operator=(static_cast<ie_loc_ghost<dim,T,layout_base,Memory>>(*this));
        (static_cast<nn_prcs<dim,T,layout_base2,Memory2>*>(&cart))->operator=(static_cast<nn_prcs<dim,T,layout_base,Memory>>(*this));
        ie_ghost<dim,T,Memory,layout_base> * ptr = static_cast<ie_ghost<dim,T,Memory,layout_base> *>((CartDecomposition<dim,T,Memory,layout_base,Distribution> *)this);
        (static_cast<ie_ghost<dim,T,Memory2,layout_base2>*>(&cart))->operator=(ptr->template duplicate<Memory2,layout_base2>());
 
        cart.private_get_sub_domains() = sub_domains;
        cart.private_get_box_nn_processor() = box_nn_processor;
        cart.private_get_fine_s() = fine_s;
        cart.private_get_gr() = gr;
        cart.private_get_gr_dist() = gr_dist;
        cart.private_get_dist() = dist;
        cart.private_get_commCostSet() = commCostSet;
        cart.private_get_cd() = cd;
        cart.private_get_domain() = domain;
        cart.private_get_sub_domains_global() = sub_domains_global;
        for (size_t i = 0 ; i < dim ; i++)
        {cart.private_get_spacing(i) = spacing[i];};
 
        cart.private_get_ghost() = ghost;
 
        cart.private_get_bbox() = bbox;
 
        for (size_t i = 0 ; i < dim ; i++)
        {cart.private_get_bc(i) = this->bc[i];}
 
        return cart;
    }
 
    CartDecomposition<dim,T,Memory, layout_base, Distribution> & operator=(const CartDecomposition & cart)
    {
        static_cast<ie_loc_ghost<dim,T,layout_base,Memory>*>(this)->operator=(static_cast<ie_loc_ghost<dim,T,layout_base,Memory>>(cart));
        static_cast<nn_prcs<dim,T,layout_base,Memory>*>(this)->operator=(static_cast<nn_prcs<dim,T,layout_base,Memory>>(cart));
        static_cast<ie_ghost<dim,T,Memory,layout_base>*>(this)->operator=(static_cast<ie_ghost<dim,T,Memory,layout_base>>(cart));
 
        sub_domains = cart.sub_domains;
        box_nn_processor = cart.box_nn_processor;
        fine_s = cart.fine_s;
        gr = cart.gr;
        gr_dist = cart.gr_dist;
        dist = cart.dist;
        commCostSet = cart.commCostSet;
        cd = cart.cd;
        domain = cart.domain;
        sub_domains_global = cart.sub_domains_global;
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            spacing[i] = cart.spacing[i];
            magn[i] = cart.magn[i];
        };
 
        ghost = cart.ghost;
 
        bbox = cart.bbox;
 
        for (size_t i = 0 ; i < dim ; i++)
            bc[i] = cart.bc[i];
 
        return *this;
    }
 
    template<typename CartDecomposition2>
    CartDecomposition<dim,T,Memory, layout_base, Distribution> & operator=(const CartDecomposition2 & cart)
    {
        static_cast<ie_loc_ghost<dim,T,layout_base,Memory>*>(this)->operator=(static_cast<typename CartDecomposition2::ie_loc_ghost_type>(cart));
        static_cast<nn_prcs<dim,T,layout_base,Memory>*>(this)->operator=(static_cast<typename CartDecomposition2::nn_prcs_type>(cart));
        static_cast<ie_ghost<dim,T,Memory,layout_base>*>(this)->operator=(static_cast<typename CartDecomposition2::ie_ghost_type>(cart));
 
        sub_domains = cart.private_get_sub_domains();
        box_nn_processor = cart.private_get_box_nn_processor();
        fine_s = cart.private_get_fine_s();
        gr = cart.private_get_gr();
        gr_dist = cart.private_get_gr_dist();
        dist = cart.private_get_dist();
        commCostSet = cart.private_get_commCostSet();
        cd = cart.private_get_cd();
        domain = cart.private_get_domain();
        sub_domains_global = cart.private_get_sub_domains_global();
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            spacing[i] = cart.private_get_spacing(i);
            magn[i] = cart.private_get_magn(i);
        };
 
        ghost = cart.private_get_ghost();
 
        bbox = cart.private_get_bbox();
 
        for (size_t i = 0 ; i < dim ; i++)
        {bc[i] = cart.private_get_bc(i);}
 
        return *this;
    }
 
    CartDecomposition<dim,T,Memory,layout_base, Distribution> & operator=(CartDecomposition && cart)
    {
        static_cast<ie_loc_ghost<dim,T,layout_base,Memory>*>(this)->operator=(static_cast<ie_loc_ghost<dim,T,layout_base,Memory>>(cart));
        static_cast<nn_prcs<dim,T,layout_base,Memory>*>(this)->operator=(static_cast<nn_prcs<dim,T,layout_base,Memory>>(cart));
        static_cast<ie_ghost<dim,T,Memory,layout_base>*>(this)->operator=(static_cast<ie_ghost<dim,T,Memory,layout_base>>(cart));
 
        sub_domains.swap(cart.sub_domains);
        box_nn_processor.swap(cart.box_nn_processor);
        fine_s.swap(cart.fine_s);
        gr = cart.gr;
        gr_dist = cart.gr_dist;
        dist = cart.dist;
        commCostSet = cart.commCostSet;
        cd = cart.cd;
        gr_dist = cart.gr_dist;
        dist = cart.dist;
 
        domain = cart.domain;
        sub_domains_global.swap(cart.sub_domains_global);
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            spacing[i] = cart.spacing[i];
            magn[i] = cart.magn[i];
        };
 
        ghost = cart.ghost;
 
        bbox = cart.bbox;
 
        for (size_t i = 0 ; i < dim ; i++)
            bc[i] = cart.bc[i];
 
        return *this;
    }
 
    static size_t getDefaultGrid(size_t n_sub)
    {
        // Calculate the number of sub-sub-domain on
        // each dimension
        return openfpm::math::round_big_2(pow(n_sub, 1.0 / dim));
    }
 
    template<typename Mem> size_t inline processorID(const encapc<1, Point<dim,T>, Mem> & p) const
    {
        return processorID_impl(p,fine_s,sub_domains_global,getDomain(),bc);
    }
 
    size_t inline processorID(const Point<dim,T> &p) const
    {
        return processorID_impl(p,fine_s,sub_domains_global,getDomain(),bc);
    }
 
    size_t inline processorID(const T (&p)[dim]) const
    {
        return processorID_impl(p,fine_s,sub_domains_global,getDomain(),bc);
    }
 
    template<typename Mem> size_t inline processorIDBC(encapc<1, Point<dim,T>, Mem> p)
    {
        Point<dim,T> pt = p;
        applyPointBC(pt);
 
 
        return processorID_impl(pt,fine_s,sub_domains_global,getDomain(),bc);
    }
 
    size_t inline processorIDBC(const Point<dim,T> &p) const
    {
        Point<dim,T> pt = p;
        applyPointBC(pt);
 
        // Get the number of elements in the cell
 
        return processorID_impl(pt,fine_s,sub_domains_global,getDomain(),bc);
    }
 
    size_t inline processorIDBC(const T (&p)[dim]) const
    {
        Point<dim,T> pt = p;
        applyPointBC(pt);
 
        return processorID_impl(pt,fine_s,sub_domains_global,getDomain(),bc);
    }
 
    inline size_t periodicity(size_t i) const
    {
        return bc[i];
    }
 
    inline const size_t (& periodicity() const) [dim]
    {
        return bc;
    }
 
    void calculate_magn(const grid_sm<dim,void> & gm)
    {
        if (gm.size() == 0)
        {
            for (size_t i = 0 ; i < dim ; i++)
                magn[i] = 1;
        }
        else
        {
            for (size_t i = 0 ; i < dim ; i++)
            {
                if (gr.size(i) % gm.size(i) != 0)
                    std::cerr << __FILE__ << ":" << __LINE__ << ".Error the decomposition grid specified as gr.size(" << i << ")=" << gr.size(i) << " is not multiple of the distribution grid gm.size(" << i << ")=" << gm.size(i) << std::endl;
 
                magn[i] = gr.size(i) / gm.size(i);
            }
        }
    }
 
    void setGoodParameters(Box<dim,T> & domain_,
                           const size_t (& bc)[dim],
                           const Ghost<dim,T> & ghost,
                           size_t dec_gran,
                           const grid_sm<dim,void> & sec_dist = grid_sm<dim,void>())
    {
        size_t div[dim];
 
        // Create a valid decomposition of the space
        // Get the number of processor and calculate the number of sub-domain
        // for decomposition
        size_t n_proc = v_cl.getProcessingUnits();
        size_t n_sub = n_proc * dec_gran;
 
        // Calculate the maximum number (before merging) of sub-domain on
        // each dimension
 
        nsub_to_div2(div,n_sub,dim);
 
/*      for (size_t i = 0; i < dim; i++)
        {
            div[i] = openfpm::math::round_big_2(pow(n_sub, 1.0 / dim));
        }*/
 
        if (dim > 3)
        {
            long int dim_r = dim-1;
            do
            {
                // Check for adjustment
                size_t tot_size = 1;
                for (size_t i = 0 ; i < dim ; i++)
                {tot_size *= div[i];}
 
                // the granularity is too coarse increase the divisions
                if (tot_size / n_proc > (unsigned int long)(0.75*dec_gran) )
                {break;}
 
                nsub_to_div(div,n_sub,dim_r);
 
                dim_r--;
            } while(dim_r > 0);
        }
 
        setParameters(div,domain_,bc,ghost,sec_dist);
    }
 
    void getParameters(size_t (& div_)[dim])
    {
        for (size_t i = 0 ; i < dim ; i++)
        {div_[i] = this->gr.size(i);}
    }
 
    void setParameters(const size_t (& div_)[dim],
                       Box<dim,T> & domain_,
                        const size_t (& bc)[dim],
                        const Ghost<dim,T> & ghost,
                        const grid_sm<dim,void> & sec_dist = grid_sm<dim,void>())
    {
        // set the boundary conditions
        for (size_t i = 0 ; i < dim ; i++)
            this->bc[i] = bc[i];
 
        // set the ghost
        this->ghost = ghost;
 
        // Set the decomposition parameters
        gr.setDimensions(div_);
        domain = domain_;
        cd.setDimensions(domain, div_, 0);
 
        // We we have a secondary grid costruct a reduced graph
        if (sec_dist.size(0) != 0)
        {
            calculate_magn(sec_dist);
            gr_dist.setDimensions(sec_dist.getSize());
        }
        else
        {
            calculate_magn(sec_dist);
            gr_dist = gr;
        }
 
        // init distribution
        dist.createCartGraph(gr_dist, domain);
 
    }
 
    void reset()
    {
        sub_domains.clear();
        box_nn_processor.clear();
        fine_s.clear();
        loc_boxes.clear();
        nn_prcs<dim, T,layout_base,Memory>::reset();
        ie_ghost<dim,T,Memory,layout_base>::reset();
        ie_loc_ghost<dim, T,layout_base,Memory>::reset();
    }
 
    void decompose(dec_options opt = dec_options::DEC_NONE)
    {
        reset();
 
        if (commCostSet == false)
        {computeCommunicationAndMigrationCosts(1);}
 
        dist.decompose();
 
        createSubdomains(v_cl,bc);
 
        calculateGhostBoxes();
 
        domain_nn_calculator_cart<dim>::reset();
        domain_nn_calculator_cart<dim>::setParameters(proc_box);
 
        if (opt != dec_options::DEC_SKIP_ICELL)
        {
 
            domain_icell_calculator<dim,T,layout_base,Memory>
			::CalculateInternalCells(v_cl,
                                 ie_ghost<dim, T,Memory,layout_base>::private_get_vb_int_box(),
                                 sub_domains,
                                 this->getProcessorBounds(),
                                 this->getGhost().getRcut(),
                                 this->getGhost());
        }
    }
 
    void refine(size_t ts)
    {
        reset();
 
        if (commCostSet == false)
        {computeCommunicationAndMigrationCosts(ts);}
 
        dist.refine();
 
        createSubdomains(v_cl,bc);
 
        calculateGhostBoxes();
 
        domain_nn_calculator_cart<dim>::reset();
        domain_nn_calculator_cart<dim>::setParameters(proc_box);
    }
 
    void redecompose(size_t ts)
    {
        reset();
 
        if (commCostSet == false)
        {computeCommunicationAndMigrationCosts(ts);}
 
        dist.redecompose();
 
        createSubdomains(v_cl,bc);
 
        calculateGhostBoxes();
 
        domain_nn_calculator_cart<dim>::reset();
        domain_nn_calculator_cart<dim>::setParameters(proc_box);
    }
 
    bool refine(DLB & dlb)
    {
        // if the DLB heuristic to use is the "Unbalance Threshold" get unbalance percentage
        if (dlb.getHeurisitc() == DLB::Heuristic::UNBALANCE_THRLD)
        {
            float unbalance = dist.getUnbalance();
            dlb.setUnbalance(unbalance);
            if (v_cl.getProcessUnitID() == 0)
            {
                std::cout << std::setprecision(3) << unbalance << "\n";
            }
        }
 
        if (dlb.rebalanceNeeded())
        {
            refine(dlb.getNTimeStepSinceDLB());
 
            return true;
        }
        return false;
    }
 
//  size_t n_step = 0;
 
    float getUnbalance()
    {
        return dist.getUnbalance();
    }
 
    size_t getProcessorLoad()
    {
        return dist.getProcessorLoad();
    }
 
    inline void getSubSubDomainPosition(size_t id, T (&pos)[dim])
    {
        dist.getSubSubDomainPosition(id, pos);
    }
 
    //TODO fix in Parmetis distribution to get only the right amount of vertices
    size_t getNSubSubDomains()
    {
        return dist.getNSubSubDomains();
    }
 
    inline void setSubSubDomainComputationCost(size_t id, size_t weight)
    {
        dist.setComputationCost(id, weight);
    }
 
    inline size_t getSubSubDomainComputationCost(size_t id)
    {
        return dist.getSubSubDomainComputationCost(id);
    }
 
    size_t subSize()
    {
        return dist.subSize();
    }
 
    size_t getNSubDomain()
    {
        return sub_domains.size();
    }
 
    Box<dim, T> getSubDomain(size_t lc)
    {
        // Create a space box
        Box<dim, T> sp;
 
        // fill the space box
 
        for (size_t k = 0; k < dim; k++)
        {
            // create the Box Low and High
            sp.setLow(k, sub_domains.template get<Box<dim,T>::p1>(lc)[k]);
            sp.setHigh(k, sub_domains.template get<Box<dim,T>::p2>(lc)[k]);
        }
 
        return sp;
    }
 
    Box<dim, T> getSubDomainWithGhost(size_t lc)
    {
        // Create a space box
        Box<dim, T> sp = sub_domains.get(lc);
 
        // enlarge with ghost
        sp.enlarge(ghost);
 
        return sp;
    }
 
    const Box<dim,T> & getDomain() const
    {
        return domain;
    }
 
    const openfpm::vector<Box<dim, T>,Memory,layout_base> &
    getSubDomains() const
    {
        return sub_domains;
    }
 
    template<typename Mem> bool isLocal(const encapc<1, Point<dim, T>, Mem> p) const
    {
        return processorID<Mem>(p) == v_cl.getProcessUnitID();
    }
 
    bool isLocal(const T (&pos)[dim]) const
    {
        return processorID(pos) == v_cl.getProcessUnitID();
    }
 
    bool isLocal(const Point<dim,T> & pos) const
    {
        return processorID(pos) == v_cl.getProcessUnitID();
    }
 
    template<typename Mem> bool isLocalBC(const encapc<1, Point<dim,T>, Mem> p, const size_t (& bc)[dim]) const
    {
        Point<dim,T> pt = p;
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (bc[i] == PERIODIC)
                pt.get(i) = openfpm::math::periodic_l(p.template get<0>()[i],domain.getHigh(i),domain.getLow(i));
        }
 
        return processorID<Mem>(pt) == v_cl.getProcessUnitID();
    }
 
    bool isLocalBC(const Point<dim,T> & p, const size_t (& bc)[dim]) const
    {
        Point<dim,T> pt = p;
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (bc[i] == PERIODIC)
                pt.get(i) = openfpm::math::periodic_l(p[i],domain.getHigh(i),domain.getLow(i));
        }
 
        return processorID(pt) == v_cl.getProcessUnitID();
    }
 
    openfpm::vector<size_t> & getDomainCells()
    {
        return domain_nn_calculator_cart<dim>::getDomainCells();
    }
 
    openfpm::vector<size_t> & getCRSDomainCells()
    {
        return domain_nn_calculator_cart<dim>::getCRSDomainCells();
    }
 
    void setNNParameters(grid_key_dx<dim> & shift, grid_sm<dim,void> & gs)
    {
        domain_nn_calculator_cart<dim>::setNNParameters(loc_boxes, shift, gs);
    }
 
    openfpm::vector<subsub_lin<dim>> & getCRSAnomDomainCells()
    {
        return domain_nn_calculator_cart<dim>::getCRSAnomDomainCells();
    }
 
    bool isLocalBC(const T (&p)[dim], const size_t (& bc)[dim]) const
    {
        Point<dim,T> pt = p;
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (bc[i] == PERIODIC)
                pt.get(i) = openfpm::math::periodic_l(p[i],domain.getHigh(i),domain.getLow(i));
        }
 
        return processorID(pt) == v_cl.getProcessUnitID();
    }
 
 
    Box<dim, T> & getProcessorBounds()
    {
        return bbox;
    }
 
 
    const Ghost<dim,T> & getGhost() const
    {
        return ghost;
    }
 
    const grid_sm<dim,void> getGrid()
    {
        return gr;
    }
 
    const grid_sm<dim,void> getDistGrid()
    {
        return gr_dist;
    }
 
 
    bool write(std::string output) const
    {
        VTKWriter<openfpm::vector<Box<dim, T>,Memory,layout_base>, VECTOR_BOX> vtk_box1;
        vtk_box1.add(sub_domains);
        vtk_box1.write(output + std::string("subdomains_") + std::to_string(v_cl.getProcessUnitID()) + std::string(".vtk"));
 
        nn_prcs<dim, T,layout_base,Memory>::write(output);
        ie_ghost<dim,T,Memory,layout_base>::write(output, v_cl.getProcessUnitID());
        ie_loc_ghost<dim, T,layout_base,Memory>::write(output, v_cl.getProcessUnitID());
 
        return true;
    }
 
    Vcluster<> & getVC() const
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif
        return v_cl;
    }
 
    void free_geo_cell()
    {
        ie_ghost<dim,T,Memory,layout_base>::free_geo_cell();
    }
 
    void free_fines()
    {
        fine_s.clear();
        fine_s.destroy();
    }
 
    bool check_consistency()
    {
        if (ie_loc_ghost<dim, T,layout_base,Memory>::check_consistency(getNSubDomain()) == false)
            return false;
 
        return true;
    }
 
    void debugPrint()
    {
        std::cout << "Subdomains\n";
        for (size_t p = 0; p < sub_domains.size(); p++)
        {
            std::cout << Box<dim, T>(sub_domains.get(p)).toString() << "\n";
        }
 
        std::cout << "Subdomains global\n";
        for (size_t p = 0; p < sub_domains_global.size(); p++)
        {
            std::cout << Box<dim, T>(sub_domains_global.template get<0>(p)).toString() << " proc:" << sub_domains_global.template get<1>(p) << "\n";
        }
 
        std::cout << "External ghost box\n";
 
        for (size_t p = 0; p<nn_prcs < dim, T, layout_base, Memory>::getNNProcessors(); p++)
        {
            for (size_t i = 0; i<ie_ghost <dim,T,Memory,layout_base>::getProcessorNEGhost(p); i++)
            {
                std::cout << ie_ghost<dim,T,Memory,layout_base>::getProcessorEGhostBox(p, i).toString() << "   prc=" << nn_prcs<dim, T, layout_base,Memory>::IDtoProc(p) << "   id=" << ie_ghost<dim,T,Memory,layout_base>::getProcessorEGhostId(p, i) << "\n";
            }
        }
 
        std::cout << "Internal ghost box\n";
 
        for (size_t p = 0; p<nn_prcs < dim, T, layout_base, Memory>::getNNProcessors(); p++)
        {
            for (size_t i = 0; i<ie_ghost<dim,T,Memory,layout_base>::getProcessorNIGhost(p); i++)
            {
                std::cout << ie_ghost<dim,T,Memory,layout_base>::getProcessorIGhostBox(p, i).toString() << "   prc=" << nn_prcs<dim, T, layout_base, Memory>::IDtoProc(p) << "   id=" << ie_ghost<dim,T,Memory,layout_base>::getProcessorIGhostId(p, i) << "\n";
            }
        }
    }
 
    bool is_equal(CartDecomposition<dim,T,Memory> & cart)
    {
        if (static_cast<ie_loc_ghost<dim,T,layout_base,Memory>*>(this)->is_equal(static_cast<ie_loc_ghost<dim,T,layout_base,Memory>&>(cart)) == false)
            return false;
 
        if (static_cast<nn_prcs<dim,T,layout_base,Memory>*>(this)->is_equal(static_cast<nn_prcs<dim,T, layout_base,Memory>&>(cart)) == false)
            return false;
 
        if (static_cast<ie_ghost<dim,T,Memory,layout_base>*>(this)->is_equal(static_cast<ie_ghost<dim,T,Memory,layout_base>&>(cart)) == false)
            return false;
 
        if (sub_domains != cart.sub_domains)
            return false;
 
        if (box_nn_processor != cart.box_nn_processor)
            return false;
 
        if (fine_s != cart.fine_s)
            return false;
 
        if (gr != cart.gr)
            return false;
 
        if (cd != cart.cd)
            return false;
 
        if (domain != cart.domain)
            return false;
 
        if (meta_compare<T[dim]>::meta_compare_f(cart.spacing,spacing) == false)
            return false;
 
        if (ghost != cart.ghost)
            return false;
 
        return true;
    }
 
    bool is_equal_ng(CartDecomposition<dim,T,Memory> & cart)
    {
        if (static_cast<ie_loc_ghost<dim,T,layout_base,Memory>*>(this)->is_equal_ng(static_cast<ie_loc_ghost<dim,T,layout_base,Memory>&>(cart)) == false)
            return false;
 
        if (static_cast<nn_prcs<dim,T,layout_base,Memory>*>(this)->is_equal(static_cast<nn_prcs<dim,T,layout_base,Memory>&>(cart)) == false)
            return false;
 
        if (static_cast<ie_ghost<dim,T,Memory,layout_base>*>(this)->is_equal_ng(static_cast<ie_ghost<dim,T,Memory,layout_base>&>(cart)) == false)
            return false;
 
        if (sub_domains != cart.sub_domains)
            return false;
 
        if (box_nn_processor != cart.box_nn_processor)
            return false;
 
        if (fine_s != cart.fine_s)
            return false;
 
        if (gr != cart.gr)
            return false;
 
        if (cd != cart.cd)
            return false;
 
        if (domain != cart.domain)
            return false;
 
        if (meta_compare<T[dim]>::meta_compare_f(cart.spacing,spacing) == false)
            return false;
 
        return true;
    }
 
    Distribution & getDistribution()
    {
        return dist;
    }
 
 
    inline void addComputationCost(size_t gid, size_t i)
    {
        size_t c = dist.getSubSubDomainComputationCost(gid);
 
        dist.setComputationCost(gid, c + i);
    }
 
    size_t get_ndec()
    {
        return dist.get_ndec();
    }
 
    const CellDecomposer_sm<dim, T, shift<dim,T>> & getCellDecomposer()
    {
        return cd;
    }
 
 
    CartDecomposition_gpu<dim,T,Memory,layout_base> toKernel()
    {
        if (host_dev_transfer == false)
        {
            fine_s.hostToDevice();
            sub_domains_global.template hostToDevice<0,1>();
            host_dev_transfer = true;
        }
 
        int bc_[dim];
 
        for (int i = 0 ; i < dim ; i++) {bc_[i] = this->periodicity(i);}
 
        CartDecomposition_gpu<dim,T,Memory,layout_base> cdg(fine_s.toKernel(),
                                                            ie_ghost<dim,T,Memory,layout_base>::toKernel(),
                                                sub_domains_global.toKernel(),
                                                getDomain(),
                                                bc_);
 
        return cdg;
    }
 
    friend extended_type;
 
    openfpm::vector<Box<dim, T>> & private_get_sub_domains()
    {
        return sub_domains;
    }
 
    const openfpm::vector<Box<dim, T>> & private_get_sub_domains() const
    {
        return sub_domains;
    }
 
    openfpm::vector<openfpm::vector<long unsigned int> > & private_get_box_nn_processor()
    {
        return box_nn_processor;
    }
 
    const openfpm::vector<openfpm::vector<long unsigned int> > & private_get_box_nn_processor() const
    {
        return box_nn_processor;
    }
 
    CellList<dim,T,Mem_fast<Memory,int>,shift<dim,T>> & private_get_fine_s()
    {
        return fine_s;
    }
 
    const CellList<dim,T,Mem_fast<Memory,int>,shift<dim,T>> & private_get_fine_s() const
    {
        return fine_s;
    }
 
    grid_sm<dim, void> & private_get_gr()
    {
        return gr;
    }
 
    const grid_sm<dim, void> & private_get_gr() const
    {
        return gr;
    }
 
    grid_sm<dim, void> & private_get_gr_dist()
    {
        return gr_dist;
    }
 
    const grid_sm<dim, void> & private_get_gr_dist() const
    {
        return gr_dist;
    }
 
    Distribution & private_get_dist()
    {
        return dist;
    }
 
    const Distribution & private_get_dist() const
    {
        return dist;
    }
 
    bool & private_get_commCostSet()
    {
        return commCostSet;
    }
 
    const bool & private_get_commCostSet() const
    {
        return commCostSet;
    }
 
    CellDecomposer_sm<dim, T, shift<dim,T>> & private_get_cd()
    {
        return cd;
    }
 
    const CellDecomposer_sm<dim, T, shift<dim,T>> & private_get_cd() const
    {
        return cd;
    }
 
    Box<dim,T> & private_get_domain()
    {
        return domain;
    }
 
    const Box<dim,T> & private_get_domain() const
    {
        return domain;
    }
 
    openfpm::vector<Box_map<dim, T>,Memory,layout_base> & private_get_sub_domains_global()
    {
        return sub_domains_global;
    }
 
    const openfpm::vector<Box_map<dim, T>,Memory,layout_base> & private_get_sub_domains_global() const
    {
        return sub_domains_global;
    }
 
    T & private_get_spacing(int i)
    {
        return spacing[i];
    }
 
    const T & private_get_spacing(int i) const
    {
        return spacing[i];
    }
 
    T & private_get_magn(int i)
    {
        return spacing[i];
    }
 
    const T & private_get_magn(int i) const
    {
        return spacing[i];
    }
 
 
    Ghost<dim,T> & private_get_ghost()
    {
        return ghost;
    }
 
    const Ghost<dim,T> & private_get_ghost() const
    {
        return ghost;
    }
 
    Box<dim,T> & private_get_bbox()
    {
        return bbox;
    }
 
    const Box<dim,T> & private_get_bbox() const
    {
        return bbox;
    }
 
    size_t & private_get_bc(int i)
    {
        return bc[i];
    }
 
    const size_t & private_get_bc(int i) const
    {
        return bc[i];
    }
};
 
 
#endif