SpaceDistributionWeight.hpp

/*
 * SpaceDistributionWeight.hpp
 *
 *  Created on: Apr 5, 2017
 *      Author: i-bird
 */

#ifndef SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTIONWEIGHT_HPP_
#define SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTIONWEIGHT_HPP_

#if 0

#include "util/mathutil.hpp"
#include "NN/CellList/CellDecomposer.hpp"

template<unsigned int dim, typename T>
class SpaceDistributionWeight
{
    Vcluster & v_cl;

    grid_sm<dim, void> gr;

    Box<dim, T> domain;

    Graph_CSR<nm_v, nm_e> gp;


    openfpm::vector<unsigned int> dec;

public:

    SpaceDistributionWeight(Vcluster & v_cl)
    :v_cl(v_cl)
    {
    }

    SpaceDistributionWeight(SpaceDistributionWeight<dim,T> && pm)
    {
        this->operator=(pm);
    }

    void createCartGraph(grid_sm<dim, void> & grid, Box<dim, T> dom)
    {
        size_t bc[dim];

        for (size_t i = 0 ; i < dim ; i++)
            bc[i] = NON_PERIODIC;

        // Set grid and domain
        gr = grid;
        domain = dom;

        // Create a cartesian grid graph
        CartesianGraphFactory<dim, Graph_CSR<nm_v, nm_e>> g_factory_part;
        gp = g_factory_part.template construct<NO_EDGE, nm_v::id, T, dim - 1, 0>(gr.getSize(), domain, bc);

        // Init to 0.0 axis z (to fix in graphFactory)
        if (dim < 3)
        {
            for (size_t i = 0; i < gp.getNVertex(); i++)
                gp.vertex(i).template get<nm_v::x>()[2] = 0.0;
        }
        for (size_t i = 0; i < gp.getNVertex(); i++)
            gp.vertex(i).template get<nm_v::global_id>() = i;

    }

    Graph_CSR<nm_v, nm_e> & getGraph()
    {
        return gp;
    }

    void decompose()
    {
        // collect all the weight on the graph index + weight on master

        openfpm::vector<id_w> id_w;
        v_cl.SGather(id_w,0);

        // Calculate the total weight
        v_cl.sum();
        v_cl.execute();

        tot /= v_cl.getProcessingUnits();

        // only master compute
        if (v_cl.getProcessUnitID() == 0)
        {
            for (size_t i = 0 ; i < id_w.size() ; i++)
                gp.template vertex_p<nm_v::computation>(id_w.get(i).id) = id_w.get(i).w;

            // Get the maximum along dimensions and take the smallest n number
            // such that 2^n < m. n it will be order of the hilbert curve

            size_t max = 0;

            for (size_t i = 0; i < dim ; i++)
            {
                if (max < gr.size(i))
                    max = gr.size(i);
            }

            // Get the order of the hilbert-curve
            size_t order = openfpm::math::log2_64(max);
            if (1ul << order < max)
                order += 1;

            size_t n = 1 << order;

            // Create the CellDecomoser

            CellDecomposer_sm<dim,T> cd_sm;
            cd_sm.setDimensions(domain, gr.getSize(), 0);

            // create the hilbert curve

            //hilbert curve iterator
            grid_key_dx_iterator_hilbert<dim> h_it(order);

            T spacing[dim];

            // Calculate the hilbert curve spacing
            for (size_t i = 0 ; i < dim ; i++)
                spacing[i] = (domain.getHigh(i) - domain.getLow(i)) / n;

            // Small grid to detect already assigned sub-sub-domains
            grid_cpu<dim,aggregate<long int>> g(gr.getSize());
            g.setMemory();

            // Reset the grid to -1
            grid_key_dx_iterator<dim> it(gr);
            while (it.isNext())
            {
                auto key = it.get();

                g.template get<0>(key) = -1;

                ++it;
            }

            // Go along the hilbert-curve and divide the space
            size_t proc_d = 0;
            size_t ele_d = 0;
            while (h_it.isNext())
            {
              auto key = h_it.get();

              // Point p
              Point<dim,T> p;

              for (size_t i = 0 ; i < dim ; i++)
                  p.get(i) = key.get(i) * spacing[i] + spacing[i] / 2;

              grid_key_dx<dim> sp = cd_sm.getCellGrid(p);

              if (g.template get<0>(sp) == -1)
              {
                  g.template get<0>(sp) = proc_d;
                  ele_d++;

                  if (ele_d >= avg)
                      proc_d++;
              }

              ++h_it;
            }

            // Fill from the grid to the graph

            // Reset the grid to -1
            grid_key_dx_iterator<dim> it2(gr);
            while (it2.isNext())
            {
                auto key = it2.get();

                gp.template vertex_p<nm_v::proc_id>(gr.LinId(key)) = g.template get<0>(key);

                id_w.get(g.template get<1>(key)) = g.template get<1>(key);

                ++it2;
            }
        }

        // Broadcast the result to all the processors.

        v_cl.Bcast(id_w);
        v_cl.execute();

        //

        return;
    }

    void refine()
    {
        decompose();
    }

    float getUnbalance()
    {
        return gr.size() % v_cl.getProcessingUnits();
    }

    void getSubSubDomainPosition(size_t id, T (&pos)[dim])
    {
#ifdef SE_CLASS1
        if (id >= gp.getNVertex())
            std::cerr << __FILE__ << ":" << __LINE__ << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex() << ")\n";
#endif

        // Copy the geometrical informations inside the pos vector
        pos[0] = gp.vertex(id).template get<nm_v::x>()[0];
        pos[1] = gp.vertex(id).template get<nm_v::x>()[1];
        if (dim == 3)
            pos[2] = gp.vertex(id).template get<nm_v::x>()[2];
    }

    inline void setComputationCost(size_t id, size_t weight)
    {
        std::cout << __FILE__ << ":" << __LINE__ << " You are trying to set the computation cost on a fixed decomposition, this operation has no effect" << std::endl;
    }

    bool weightsAreUsed()
    {
        return false;
    }

    size_t getSubSubDomainComputationCost(size_t id)
    {
        return 1.0;
    }

    size_t getProcessorLoad()
    {
        // Get the number of processing units
        size_t Np = v_cl.getProcessingUnits();

        // Calculate the best number of sub-domains for each
        // processor
        size_t N_tot = gr.size();
        size_t N_best_each = N_tot / Np;
        size_t N_rest = N_tot % Np;

        if (v_cl.getProcessUnitID() < N_rest)
            N_best_each += 1;

        return N_best_each;
    }

    void setMigrationCost(size_t id, size_t migration)
    {
    }

    void setCommunicationCost(size_t v_id, size_t e, size_t communication)
    {
    }

    size_t getNSubSubDomains()
    {
        return gp.getNVertex();
    }

    size_t getNSubSubDomainNeighbors(size_t id)
    {
        return gp.getNChilds(id);
    }

    void write(const std::string & file)
    {
        VTKWriter<Graph_CSR<nm_v, nm_e>, VTK_GRAPH> gv2(gp);
        gv2.write(std::to_string(v_cl.getProcessUnitID()) + "_" + file + ".vtk");
    }

    const SpaceDistribution<dim,T> & operator=(const SpaceDistribution<dim,T> & dist)
    {
        gr = dist.gr;
        domain = dist.domain;
        gp = dist.gp;

        return *this;
    }

    const SpaceDistribution<dim,T> & operator=(SpaceDistribution<dim,T> && dist)
    {
        v_cl = dist.v_cl;
        gr = dist.gr;
        domain = dist.domain;
        gp.swap(dist.gp);

        return *this;
    }

    size_t get_ndec()
    {
        return 0;
    }
};

#endif

#endif /* SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTIONWEIGHT_HPP_ */