HyperCube.hpp

#ifndef HYPERCUBE_HPP
#define HYPERCUBE_HPP
 
#include "Grid/grid_sm.hpp"
#include "Grid/comb.hpp"
#include "util/mathutil.hpp"
 
template<unsigned int dim, unsigned int subdim> class SubHyperCube;
 
template<unsigned int dim> std::ostream& operator<<(std::ostream& str, const comb<dim> & c)
{
    // print the combination of ostream
    for (size_t i = 0 ; i < dim-1 ; i++)
    {
        /* coverity[dead_error_line] */
        str <<  c.c[i] << ";";
    }
 
    str << c.c[dim-1];
 
    return str;
}
 
template<unsigned int dim>
class HyperCube
{
public:
 
    static size_t getNumberOfElements_R(size_t d)
    {
        // The formula to calculate the number of element of size d is given by
        //
        // 2^(dim-d) * C(dim,d) with C(dim,d) = dim!/(d!(dim-d)!)
 
        return pow(2,dim-d) * openfpm::math::C(dim,d);
    }
 
    static size_t getNumberOfElementsTo_R(size_t d_t)
    {
        size_t N_ele = 0;
 
        for (size_t i = 0 ; i <= d_t ; i++)
            N_ele += getNumberOfElements_R(i);
 
        return N_ele;
    }
 
    static std::vector<comb<dim>> getCombinations_R(size_t d)
    {
#ifdef SE_CLASS1
        if (d > dim)
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " " << d << " must be smaller than " << "\n";
#endif
 
        // Create an Iterator_g_const
        // And a vector that store all the combination
 
        std::vector<comb<dim>> v;
        Iterator_g_const it(dim-d,dim);
 
        // for each non-zero elements
        // basically the key will store the position of the
        // non zero elements, while BinPermutations will
        // fill the array of all the permutations
        //
 
        while (it.isNext())
        {
            grid_key_dx_r & key = it.get();
 
            // Calculate the permutation
 
            BinPermutations(key,v);
            ++it;
        }
 
        // case when d == dim
 
        if (d == dim)
        {
            comb<dim> c;
            c.zero();
 
            v.push_back(c);
        }
 
        // return the combinations
 
        return v;
    }
 
    static std::vector<comb<dim>> getCombinations_R_bc(size_t d, const size_t (& bc)[dim])
    {
#ifdef SE_CLASS1
        if (d > dim)
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " " << d << " must be smaller than " << "\n";
#endif
 
        size_t dg[dim];
 
        size_t k = 0;
        // get the indexes of the free degree of freedom
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (bc[i] == PERIODIC)
            {
                dg[k] = i;
                k++;
            }
        }
 
        // Create an Iterator_g_const
        // And a vector that store all the combination
 
        std::vector<comb<dim>> v;
        Iterator_g_const it(k-(d-(dim - k)),k);
 
        // for each non-zero elements
        // basically the key will store the position of the
        // non zero elements, while BinPermutations will
        // fill the array of all the permutations
        //
 
        while (it.isNext())
        {
            grid_key_dx_r key = it.get();
 
            // Now we adjust the non zero position
 
            for (size_t i = 0 ; i < key.getDim() ; i++)
            {key.set_d(i,dg[key.get(i)]);}
 
            // Calculate the permutation
 
            BinPermutations(key,v);
            ++it;
        }
 
        // case when d == dim
 
        if (d == dim)
        {
            comb<dim> c;
            c.zero();
 
            v.push_back(c);
        }
 
        // return the combinations
 
        return v;
    }
 
    static void BinPermutations(grid_key_dx_r & pos, std::vector<comb<dim>> & v)
    {
        size_t p_stop = pow(2,pos.getDim());
        for (size_t p = 0 ; p < p_stop ; p++)
        {
            // Create a new permutation
            struct comb<dim> c;
 
            // zero the combination
            c.zero();
 
            // the bit of p give the permutations 0 mean bottom 1 mean up
            // Fill the combination based on the bit of p
 
            for (size_t k = 0 ; k < pos.getDim() ; k++)
            {
                // bit of p
                bool bit = (p >> k) & 0x01;
 
                // Fill the combination
                c.c[pos.get(k)] = (bit == 0)? 1 : -1;
            }
 
            // save the combination
 
            v.push_back(c);
        }
    }
 
    static void BinPermutationsSt(std::vector<comb<dim>> & v)
    {
        size_t p_stop = pow(2,dim);
        for (size_t p = 0 ; p < p_stop ; p++)
        {
            // Create a new permutation
            struct comb<dim> c;
 
            // zero the combination
            c.zero();
 
            // the bit of p give the permutations 0 mean bottom 1 mean up
            // Fill the combination based on the bit of p
 
            for (size_t k = 0 ; k < dim ; k++)
            {
                // bit of p
                bool bit = (p >> k) & 0x01;
 
                // Fill the combination
                c.c[k] = (bit == 0)? 0 : -1;
            }
 
            // save the combination
 
            v.push_back(c);
        }
    }
 
    static size_t LinPerm(comb<dim> & c)
    {
        size_t id = 0;
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (c.c[i] == -1)
            {id = id | (1 << i);}
        }
 
        return id;
    }
 
    /*
    static SubHyperCube<dim,dim-1> getSubHyperCube(int d)
    {
        SubHyperCube<dim,dim-1> sub;
 
        return sub;
    }*/
 
    static size_t LinId(comb<dim> & c)
    {
        // calculate the numbers of non-zero
        size_t d = 0;
 
        size_t pos_n_zero[dim];
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (c.c[i] != 0)
            {d++;}
        }
 
        // Get the position of the non-zero
        size_t pn_zero = 0;
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (c.c[i] != 0)
            {
                pos_n_zero[d-pn_zero-1] = i;
                pn_zero++;
            }
        }
 
        size_t lin_id = 0;
 
        // Cumulative value
        long int val = 0;
        long int cum_val = 0;
        for(long int i = d - 1; i >= 0 ; i--)
        {
            // check the out-of-bound, outside is assumed to be -1 so  (- pos_n_zero[i+1] - 1) = 0
            if (i+1 < (long int)d)
            {
                /* coverity[dead_error_line] */
                val = pos_n_zero[i] - pos_n_zero[i+1] - 1;
            }
            else
            {
                /* coverty[uninit_use] */
                val = pos_n_zero[i];
            }
 
            for (long int j = 0 ; j < (long int)val; j++)
            {
                // C is not safe, check the limit
                /* coverity[dead_error_line] */
                if (((long int)dim)-cum_val-j-1 >= 0 && i > 0 && ((long int)dim)-cum_val-j >= i)
                    lin_id += openfpm::math::C(dim-cum_val-j-1,i);
                else
                    lin_id += 1;
            }
 
            cum_val += (val + 1);
        }
 
        // multiply for the permutation
        lin_id *= pow(2,d);
 
        // calculate the permutation position
        size_t id = 0;
 
        for (size_t i = 0 ; i < d ; i++)
        {
            if (c.c[pos_n_zero[i]] == -1)
            {id = id | (1 << i);}
        }
 
        // return the correct id
 
        return lin_id + id;
    }
 
    static bool isPositive(size_t d)
    {
        return (d % 2) == 0;
    }
 
    static int positiveFace(int d)
    {
        return d * 2;
    }
 
    static int negativeFace(int d)
    {
        return d * 2 + 1;
    }
};
 
template<unsigned int dim, unsigned int subdim>
class SubHyperCube : public HyperCube<subdim>
{
public:
 
    static std::vector<comb<dim>> getCombinations_R(comb<dim> c, int d)
    {
#ifdef DEBUG
        if (c.n_zero() < d)
        {
            std::cerr << "Error SubHyperCube: " << __FILE__ << " " << __LINE__ << " the hyper-cube selected must have dimensionality bigger than the dimensionality of the requested combinations, or the number of zero in c must me bigger than d" << "\n";
        }
#endif
 
        // if sub-dim == 0 return c
 
        if (subdim == 0)
        {
            std::vector<comb<dim>> vc;
            vc.push_back(c);
 
            return vc;
        }
 
        // Create an Iterator_g_const
        // And a vector that store all the combination
 
        std::vector<comb<subdim>> v = HyperCube<subdim>::getCombinations_R(d);
 
        // Create new combinations
        std::vector<comb<dim>> vc(v.size());
 
        // for each combination
        for (size_t i = 0 ; i < v.size() ; i++)
        {
            // sub j counter
            int sub_j = 0;
            // for each zero (direction spanned by the sub-hyper-cube)
            for (size_t j = 0 ; j < dim ; j++)
            {
                if (c.c[j] == 0)
                {
                    // take the combination from the sub-hyper-cube
                    vc[i].c[j] = v[i].c[sub_j];
                    sub_j++;
                }
                else
                {
                    // take the combination from the hyper-cube position
                    vc[i].c[j] = c.c[j];
                }
            }
        }
 
        // return the combinations
        return vc;
    }
};
 
#endif