CellList.hpp

/*
 * CellList.hpp
 *
 *  Created on: Mar 21, 2015
 *      Author: Pietro Incardona
 */

#ifndef CELLLIST_HPP_
#define CELLLIST_HPP_

#include "Vector/map_vector.hpp"
#include "CellDecomposer.hpp"
#include "Space/SpaceBox.hpp"
#include "util/mathutil.hpp"
#include "CellNNIterator.hpp"
#include "Space/Shape/HyperCube.hpp"
#include "CellListNNIteratorRadius.hpp"
#include <unordered_map>

#include "CellListIterator.hpp"
#include "ParticleIt_Cells.hpp"
#include "ParticleItCRS_Cells.hpp"
#include "util/common.hpp"

#include "NN/Mem_type/MemFast.hpp"
#include "NN/Mem_type/MemBalanced.hpp"
#include "NN/Mem_type/MemMemoryWise.hpp"
#include "NN/CellList/NNc_array.hpp"

template<typename key,typename val>
class wrap_unordered_map: public std::unordered_map<key,val>
{
};

#ifdef HAVE_LIBQUADMATH

#include <boost/multiprecision/float128.hpp>


template<typename val>
class wrap_unordered_map<boost::multiprecision::float128,val>
{
};

#endif

#define CELL_REALLOC 16ul

#define STARTING_NSLOT 16

template<unsigned int dim> void NNcalc_csr(openfpm::vector<std::pair<grid_key_dx<dim>,grid_key_dx<dim>>> & cNN)
{
    // Calculate the NNc_full array, it is a structure to get the neighborhood array

    // compile-time array {0,0,0,....}  {2,2,2,...} {1,1,1,...}

    typedef typename generate_array<size_t,dim, Fill_zero>::result NNzero;
    typedef typename generate_array<size_t,dim, Fill_two>::result NNtwo;

    // Generate the sub-grid iterator

    size_t div[dim];

    // Calculate the divisions

    for (size_t i = 0 ; i < dim ; i++)
        div[i] = 4;

    grid_sm<dim,void> gs(div);
    grid_key_dx_iterator_sub<dim> gr_sub3(gs,NNzero::data,NNtwo::data);

    grid_key_dx<dim> src_; // Source cell
    for (size_t i = 0; i < dim; i++)
        src_.set_d(i,1);

    size_t middle = gs.LinId(src_);

    // Calculate the symmetric crs array
    while (gr_sub3.isNext())
    {
        auto dst = gr_sub3.get();
        grid_key_dx<dim> src = src_;

        if ((long int)middle > gs.LinId(dst))
        {
            ++gr_sub3;
            continue;
        }

        // Here we adjust src and dst to be in the positive quadrant

        for (size_t i = 0 ; i < dim; i++)
        {
            if (dst.get(i) == 0)
            {
                src.set_d(i,src.get(i) + 1);
                dst.set_d(i,dst.get(i) + 1);
            }
        }

        src -= src_;
        dst -= src_;

        cNN.add(std::pair<grid_key_dx<dim>,grid_key_dx<dim>>(src,dst));

        ++gr_sub3;

    }
};

template<unsigned int dim> void NNcalc_sym(openfpm::vector<grid_key_dx<dim>> & cNN)
{
    // Calculate the NNc_full array, it is a structure to get the neighborhood array

    // compile-time array {0,0,0,....}  {2,2,2,...} {1,1,1,...}

    typedef typename generate_array<size_t,dim, Fill_zero>::result NNzero;
    typedef typename generate_array<size_t,dim, Fill_two>::result NNtwo;

    // Generate the sub-grid iterator

    size_t div[dim];

    // Calculate the divisions

    for (size_t i = 0 ; i < dim ; i++)
        div[i] = 4;

    grid_sm<dim,void> gs(div);
    grid_key_dx_iterator_sub<dim> gr_sub3(gs,NNzero::data,NNtwo::data);

    grid_key_dx<dim> src_; // Source cell
    for (size_t i = 0; i < dim; i++)
        src_.set_d(i,1);

    size_t middle = gs.LinId(src_);

    // Calculate the symmetric array
    while (gr_sub3.isNext())
    {
        auto dst = gr_sub3.get();

        if ((long int)middle > gs.LinId(dst))
        {
            ++gr_sub3;
            continue;
        }

        cNN.add(dst - src_);

        ++gr_sub3;

    }
};

template<unsigned int dim> void NNcalc_full(openfpm::vector<grid_key_dx<dim>> & cNN)
{
    // Calculate the NNc_full array, it is a structure to get the neighborhood array

    // compile-time array {0,0,0,....}  {2,2,2,...} {1,1,1,...}

    typedef typename generate_array<size_t,dim, Fill_zero>::result NNzero;
    typedef typename generate_array<size_t,dim, Fill_two>::result NNtwo;

    // Generate the sub-grid iterator

    size_t div[dim];

    // Calculate the divisions

    for (size_t i = 0 ; i < dim ; i++)
    {div[i] = 4;}

    grid_sm<dim,void> gs(div);
    grid_key_dx_iterator_sub<dim> gr_sub3(gs,NNzero::data,NNtwo::data);

    grid_key_dx<dim> src_; // Source cell
    for (size_t i = 0; i < dim; i++)
        src_.set_d(i,1);

    // Calculate the symmetric crs array
    while (gr_sub3.isNext())
    {
        auto dst = gr_sub3.get();

        cNN.add(dst - src_);

        ++gr_sub3;
    }
};


/* NOTE all the implementations
 *
 * has complexity O(1) in getting the cell id and the elements in a cell
 * but with different constants
 *
 */

template<unsigned int dim, typename T,  typename Mem_type, typename transform = no_transform<dim,T>, typename base=openfpm::vector<size_t>>
class CellList : public CellDecomposer_sm<dim,T,transform>, public Mem_type
{
protected:
    //
    //    * * *
    //    * x *
    //    * * *


    NNc_array<dim,(unsigned int)openfpm::math::pow(3,dim)> NNc_full;
//  long int NNc_full[openfpm::math::pow(3,dim)];

    //
    //   * * *
    //     x *
    //
//  long int NNc_sym[openfpm::math::pow(3,dim)/2+1];
    NNc_array<dim,(unsigned int)openfpm::math::pow(3,dim)/2+1> NNc_sym;

private:

    wrap_unordered_map<T,openfpm::vector<long int>> rcache;

    bool from_cd;

    size_t n_dec;

    void NNcalc(T r_cut, openfpm::vector<long int> & NNcell)
    {
        size_t n_cell[dim];
        size_t n_cell_mid[dim];

        Point<dim,T> spacing = this->getCellBox().getP2();
        const grid_sm<dim,void> & gs = this->getGrid();

        for (size_t i = 0 ; i < dim ; i++)
        {
            n_cell[i] = 2*(std::ceil(r_cut / spacing.get(i)))+1;
            n_cell_mid[i] = n_cell[i] / 2;
        }

        grid_sm<dim,void> gsc(n_cell);
        grid_key_dx_iterator<dim> gkdi(gsc);

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

            for (size_t i = 0 ; i < dim ; i++)
                key.set_d(i,key.get(i) - n_cell_mid[i]);

            NNcell.add(gs.LinId(key));

            ++gkdi;
        }
    }

    void InitializeStructures(const size_t (& div)[dim], size_t tot_n_cell, size_t slot=STARTING_NSLOT)
    {
        Mem_type::init_to_zero(slot,tot_n_cell);

        NNc_full.set_size(div);
        NNc_full.init_full();

        NNc_sym.set_size(div);
        NNc_sym.init_sym();
    }

    void setCellDecomposer(CellDecomposer_sm<dim,T,transform> & cd, const CellDecomposer_sm<dim,T,transform> & cd_sm, const Box<dim,T> & dom_box, size_t pad) const
    {
        size_t bc[dim];

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

        Box<dim,long int> bx = cd_sm.convertDomainSpaceIntoCellUnits(dom_box,bc);

        size_t div[dim];
        size_t div_big[dim];

        for (size_t i = 0 ; i < dim ; i++)
        {
            div[i] = bx.getHigh(i) - bx.getLow(i);
            div_big[i] = cd_sm.getDiv()[i] - 2*cd_sm.getPadding(i);
        }

        cd.setDimensions(cd_sm.getDomain(),div_big,div, pad, bx.getP1());
    }

public:

    typedef Mem_type Mem_type_type;

    typedef CellNNIteratorSym<dim,CellList<dim,T,Mem_type,transform,base>,RUNTIME,NO_CHECK> SymNNIterator;

    typedef typename base::value_type value_type;

    typedef T stype;

    const grid_sm<dim,void> & getGrid()
    {
        return CellDecomposer_sm<dim,T,transform>::getGrid();
    }

    void Initialize(CellDecomposer_sm<dim,T,transform> & cd_sm, const Box<dim,T> & dom_box,const size_t pad = 1, size_t slot=STARTING_NSLOT)
    {
        size_t bc[dim];
        for (size_t i = 0 ; i < dim ; i++)  {bc[i] = NON_PERIODIC;}

        Box<dim,long int> bx = cd_sm.convertDomainSpaceIntoCellUnits(dom_box,bc);

        setCellDecomposer(*this,cd_sm,dom_box,pad);

        size_t div_w_pad[dim];
        size_t tot_cell = 1;

        for (size_t i = 0 ; i < dim ; i++)
        {
            div_w_pad[i] = bx.getHigh(i) - bx.getLow(i) + 2*pad;
            tot_cell *= div_w_pad[i];
        }

        // here we set the cell-shift for the CellDecomposer

        InitializeStructures(div_w_pad,tot_cell);

        // Initialized from CellDecomposer
        from_cd = true;
    }

    void Initialize(const Box<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
    {
        SpaceBox<dim,T> sbox(box);

        // Initialize point transformation

        Initialize(sbox,div,pad,slot);
    }

    void Initialize(const SpaceBox<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
    {
        Matrix<dim,T> mat;

        CellDecomposer_sm<dim,T,transform>::setDimensions(box,div, mat, pad);
        Mem_type::set_slot(slot);

        // create the array that store the number of particle on each cell and se it to 0
        InitializeStructures(this->gr_cell.getSize(),this->gr_cell.size());

        from_cd = false;
    }

    CellList()

    :Mem_type(STARTING_NSLOT)
    {};

    CellList(const CellList<dim,T,Mem_type,transform,base> & cell)
    :Mem_type(STARTING_NSLOT)
    {
        this->operator=(cell);
    }

    CellList(CellList<dim,T,Mem_type,transform,base> && cell)
    :Mem_type(STARTING_NSLOT)
    {
        this->operator=(cell);
    }


    CellList(Box<dim,T> & box, const size_t (&div)[dim], Matrix<dim,T> mat, const size_t pad = 1, size_t slot=STARTING_NSLOT)
    :Mem_type(slot),CellDecomposer_sm<dim,T,transform>(box,div,mat,box.getP1(),pad)
    {
        SpaceBox<dim,T> sbox(box);
        Initialize(sbox,div,pad,slot);
    }

    CellList(Box<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
    :Mem_type(slot),n_dec(0)
    {
        SpaceBox<dim,T> sbox(box);
        Initialize(sbox,div,pad,slot);
    }

    CellList(SpaceBox<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
    :Mem_type(slot)
    {
        Initialize(box,div,pad,slot);
    }

    CellList(CellDecomposer_sm<dim,T,transform> & cd_sm, const Box<dim,T> & box, const size_t pad = 1, size_t slot=STARTING_NSLOT)
    :Mem_type(slot),n_dec(0)
    {
        Initialize(cd_sm,box,pad,slot);
    }

    ~CellList()
    {}

    CellList<dim,T,Mem_type,transform,base> & operator=(CellList<dim,T,Mem_type,transform,base> && cell)
    {
        std::copy(&cell.NNc_full[0],&cell.NNc_full[openfpm::math::pow(3,dim)],&NNc_full[0]);
        std::copy(&cell.NNc_sym[0],&cell.NNc_sym[openfpm::math::pow(3,dim)/2+1],&NNc_sym[0]);

        Mem_type::swap(static_cast<Mem_type &&>(cell));

        static_cast<CellDecomposer_sm<dim,T,transform> &>(*this).swap(cell);

        n_dec = cell.n_dec;
        from_cd = cell.from_cd;

        return *this;
    }

    CellList<dim,T,Mem_type,transform,base> & operator=(const CellList<dim,T,Mem_type,transform,base> & cell)
    {
        NNc_full = cell.NNc_full;
        NNc_sym = cell.NNc_sym;

        Mem_type::operator=(static_cast<const Mem_type &>(cell));

        static_cast<CellDecomposer_sm<dim,T,transform> &>(*this) = static_cast<const CellDecomposer_sm<dim,T,transform> &>(cell);

        n_dec = cell.n_dec;
        from_cd = cell.from_cd;

        return *this;
    }

    ParticleIt_Cells<dim,CellList<dim,T,Mem_fast<>,transform,base>> getDomainIterator(openfpm::vector<size_t> & dom_cells)
    {
        ParticleIt_Cells<dim,CellList<dim,T,Mem_fast<>,transform,base>> it(*this,dom_cells);

        return it;
    }

    inline void addCell(size_t cell_id, typename base::value_type ele)
    {
        Mem_type::addCell(cell_id,ele);
    }

    inline void add(const T (& pos)[dim], typename base::value_type ele)
    {
        // calculate the Cell id
        size_t cell_id = this->getCell(pos);

        Mem_type::add(cell_id,ele);
    }

    inline void add(const Point<dim,T> & pos, typename base::value_type ele)
    {
        // calculate the Cell id
        size_t cell_id = this->getCell(pos);

        Mem_type::add(cell_id,ele);
    }


    inline void addDom(const T (& pos)[dim], typename base::value_type ele)
    {
        // calculate the Cell id

        size_t cell_id = this->getCellDom(pos);

        // add the element to the cell

        addCell(cell_id,ele);
    }

    inline void addDom(const Point<dim,T> & pos, typename base::value_type ele)
    {
        // calculate the Cell id

        size_t cell_id = this->getCellDom(pos);

        // add the element to the cell

        addCell(cell_id,ele);
    }

    inline void addPad(const T (& pos)[dim], typename base::value_type ele)
    {
        // calculate the Cell id

        size_t cell_id = this->getCellPad(pos);

        // add the element to the cell

        addCell(cell_id,ele);
    }

    inline void addPad(const Point<dim,T> & pos, typename base::value_type ele)
    {
        // calculate the Cell id

        size_t cell_id = this->getCell(pos);

        // add the element to the cell

        addCell(cell_id,ele);
    }

    inline void remove(size_t cell, size_t ele)
    {
        Mem_type::remove(cell,ele);
    }

    inline size_t getNelements(const size_t cell_id) const
    {
        return Mem_type::getNelements(cell_id);
    }

    inline auto get(size_t cell, size_t ele) -> decltype(this->Mem_type::get(cell,ele))
    {
        return Mem_type::get(cell,ele);
    }

    inline auto get(size_t cell, size_t ele) const -> decltype(this->Mem_type::get(cell,ele))
    {
        return Mem_type::get(cell,ele);
    }

    inline void swap(CellList<dim,T,Mem_type,transform,base> & cl)
    {
        NNc_full.swap(cl.NNc_full);
        NNc_sym.swap(cl.NNc_sym);

        Mem_type::swap(static_cast<Mem_type &>(cl));

        static_cast<CellDecomposer_sm<dim,T,transform> &>(*this) = static_cast<const CellDecomposer_sm<dim,T,transform> &>(cl);
    }

    CellIterator<CellList<dim,T,Mem_type,transform,base>> getCellIterator(size_t cell)
    {
        return CellIterator<CellList<dim,T,Mem_type,transform,base>>(cell,*this);
    }

    template<unsigned int impl=NO_CHECK> inline CellNNIterator<dim,CellList<dim,T,Mem_type,transform,base>,(int)FULL,impl> getNNIterator(size_t cell)
    {
        CellNNIterator<dim,CellList<dim,T,Mem_type,transform,base>,(int)FULL,impl> cln(cell,NNc_full,*this);
        return cln;

    }

    template<unsigned int impl=NO_CHECK> inline CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform,base>,impl> getNNIteratorRadius(size_t cell, T r_cut)
    {
        openfpm::vector<long int> & NNc = rcache[r_cut];

        if (NNc.size() == 0)
            NNcalc(r_cut,NNc);

        CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform,base>,impl> cln(cell,NNc,*this);

        return cln;
    }



    template<unsigned int impl>
    inline CellNNIteratorSym<dim,CellList<dim,T,Mem_type,transform,base>,(unsigned int)SYM,impl>
    getNNIteratorSym(size_t cell, size_t p, const openfpm::vector<Point<dim,T>> & v)
    {
#ifdef SE_CLASS1
        if (from_cd == false)
        {std::cerr << __FILE__ << ":" << __LINE__ << " Warning when you try to get a symmetric neighborhood iterator, you must construct the Cell-list in a symmetric way" << std::endl;}
#endif

        CellNNIteratorSym<dim,CellList<dim,T,Mem_type,transform,base>,SYM,impl> cln(cell,p,NNc_sym,*this,v);
        return cln;
    }

    template<unsigned int impl>
    inline CellNNIteratorSymMP<dim,CellList<dim,T,Mem_type,transform,base>,(unsigned int)SYM,impl>
    getNNIteratorSymMP(size_t cell, size_t p, const openfpm::vector<Point<dim,T>> & v_p1, const openfpm::vector<Point<dim,T>> & v_p2)
    {
#ifdef SE_CLASS1
        if (from_cd == false)
            std::cerr << __FILE__ << ":" << __LINE__ << " Warning when you try to get a symmetric neighborhood iterator, you must construct the Cell-list in a symmetric way" << std::endl;
#endif

        CellNNIteratorSymMP<dim,CellList<dim,T,Mem_type,transform,base>,SYM,impl> cln(cell,p,NNc_sym,*this,v_p1,v_p2);
        return cln;
    }

    const NNc_array<dim,(unsigned int)openfpm::math::pow(3,dim)/2+1> & getNNc_sym() const
    {
        return NNc_sym;
    }

    size_t getPadding(size_t i) const
    {
        return CellDecomposer_sm<dim,T,transform>::getPadding(i);
    }

    size_t (& getPadding())[dim]
    {
        return CellDecomposer_sm<dim,T,transform>::getPadding();
    }

    void clear()
    {
        Mem_type::clear();
    }

    inline const typename Mem_type::loc_index & getStartId(typename Mem_type::loc_index cell_id) const
    {
        return Mem_type::getStartId(cell_id);
    }

    inline const typename Mem_type::loc_index & getStopId(typename Mem_type::loc_index cell_id) const
    {
        return Mem_type::getStopId(cell_id);
    }

    inline const typename Mem_type::loc_index & get_lin(const typename Mem_type::loc_index * part_id) const
    {
        return Mem_type::get_lin(part_id);
    }


    size_t g_m = 0;

    inline size_t get_gm()
    {
        return g_m;
    }

    inline void set_gm(size_t g_m)
    {
        this->g_m = g_m;
    }



        void set_ndec(size_t n_dec)
        {
            this->n_dec = n_dec;
        }

        size_t get_ndec() const
        {
            return n_dec;
        }

};

template<unsigned int dim, typename St> static inline void cl_param_calculate(Box<dim, St> & pbox, size_t (&div)[dim], St r_cut, const Ghost<dim, St> & enlarge)
{
    // calculate the parameters of the cell list

    // extend by the ghost
    pbox.enlarge(enlarge);

    // Calculate the division array and the cell box
    for (size_t i = 0; i < dim; i++)
    {
        div[i] = static_cast<size_t>((pbox.getP2().get(i) - pbox.getP1().get(i)) / r_cut);
        div[i]++;
        pbox.setHigh(i,pbox.getLow(i) + div[i]*r_cut);
    }
}

template<unsigned int dim, typename St> static
inline void cl_param_calculateSym(const Box<dim,St> & dom,
                                  CellDecomposer_sm<dim,St,shift<dim,St>> & cd_sm,
                                  Ghost<dim,St> g,
                                  St r_cut,
                                  size_t & pad)
{
    size_t div[dim];

    for (size_t i = 0 ; i < dim ; i++)
        div[i] = (dom.getHigh(i) - dom.getLow(i)) / r_cut;

    g.magnify(1.013);

    // Calculate the maximum padding
    for (size_t i = 0 ; i < dim ; i++)
    {
        size_t tmp = std::ceil(fabs(g.getLow(i)) / r_cut);
        pad = (pad > tmp)?pad:tmp;

        tmp = std::ceil(fabs(g.getHigh(i)) / r_cut);
        pad = (pad > tmp)?pad:tmp;
    }

    cd_sm.setDimensions(dom,div,pad);
}


template<unsigned int dim, typename St>
static inline void cl_param_calculateSym(const Box<dim,St> & dom,
                                         CellDecomposer_sm<dim,St,shift<dim,St>> & cd_sm,
                                         Ghost<dim,St> g,
                                         size_t & pad)
{
    size_t div[dim];

    for (size_t i = 0 ; i < dim ; i++)
        div[i] = (dom.getHigh(i) - dom.getLow(i)) / g.getHigh(i);

    g.magnify(1.013);

    pad = 1;

    cd_sm.setDimensions(dom,div,pad);
}


#endif /* CELLLIST_HPP_ */