gsCombinatorics.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsLinearAlgebra.h>
#include <gsCore/gsMath.h>

namespace gismo
{

// \ingroup Utils
// could also be in Utils, but doxygen allows only one group for free functions
inline unsigned factorial( unsigned n)
{
    static const unsigned precomputed[]= {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600};
    GISMO_ASSERT(n < util::size(precomputed), "Overflow when computing factorial.");
    return precomputed[n];
}

/*
// \brief Returns the factorial of \a n
// Remember that factorial grow too fast and only n! with n<=21 can
// be stored in a 64bit that is an unsigned long long.  \ingroup
// combinatorics
inline unsigned long long factorial( unsigned long long n)
{
static const unsigned long long precomputed[]= {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200, 1307674368000, 20922789888000, 355687428096000, 6402373705728000, 121645100408832000, 2432902008176640000 };
GISMO_ASSERT(n < util::size(precomputed), "Overflow when computing factorial.");
return precomputed[n];
}
*/

// \ingroup Utils
// Note: could also be in Utils, but doxygen allows only one group for free functions
template <typename Z>
inline Z binomial(const Z n, const Z r)
{
    GISMO_ASSERT(r>=0, "binomial coefficient (n,r) exists only for n,r positive");
    //GISMO_ASSERT(n>=r, "binomial coefficient (n,r) exists only for n,r such that n>=r");

    const Z diff = math::min( n-r, r );
    int result = 1;
    for (Z i=0;i < diff;)
    {
        result *= n-i;
        result /= ++i;
    }
    return result;
}

//\ingroup Utils
inline void binomial_into( unsigned n, gsVector<unsigned>& v)
{
    v.resize (n+1);
    v[0]=1;

    for ( unsigned i = 1; i <= n; ++i)
    {
        v[i]= 1;
        for ( unsigned j= i-1; j>0; --j)
            v[j] += v[j-1];
    }
}

// Used by binomial<n,r>()
template<int n, int r> class binomialT
{public:enum { value= binomialT<n-1,r-1>::value+binomialT<n-1,r>::value};};
template<int n> class binomialT<n,n> {public:enum { value= 1};};
template<>      class binomialT<0,0> {public:enum { value= 1};};
template<int n> class binomialT<n,0> {public:enum { value= 1};};
template <unsigned n, unsigned r>
inline unsigned binomial () {return binomialT<n,r>::value;}

template<class Vec>
inline Vec stridesOf(const Vec & sz)
{
    Vec result(sz.size());
    result[0] = 1;
    for ( index_t i=1; i != sz.size(); ++i )
        result[i] = result[i-1] * sz[i-1];
    return result;
}

template<class Vec>
void firstCombination(const unsigned n, const unsigned r, Vec & res)
{
    if (r<=n) 
        res= Vec::LinSpaced(r,0,r-1);
    else
        std::cerr << "Error: r>n combination requested. r="<< r<<", n="<< n<<"\n";
}
  
template<class Vec>
bool nextCombination (Vec & v, const unsigned n)
{
    const index_t r = v.rows() ;            
    if  (v == Vec::LinSpaced(r,n-r,n-1)) return false;
    int i = r-1;
    while (v[i] == n-r+i) --i;
    v[i] += 1;
    for (index_t j=i+1; j<r; ++j)
        v[j] = v[i] +j-i;
    return true;
}

template<class Vec>
void firstPermutation (Vec &current)
{
    const index_t n=current.size();
    current=Vec::LinSpaced(n,0,n-1);
}

template<class Vec>
bool nextPermutation (Vec &current)
{
    const index_t n=current.size();
    return std::next_permutation(current.data(), current.data()+n);
}


template<class Vec>
bool nextLexicographic(Vec& cur, const Vec& size)
{
    const index_t d = cur.size();
    GISMO_ASSERT( d == size.size(), "Vector sizes don't match in nextLexicographic" );

    for (index_t i = 0; i < d; ++i)
    {
        // increase current dimension
        if (++cur[i] == size[i])    // current dimension exhausted ?
        {
            if (i == d - 1)         // was it the last one?
                return false;       // then all elements exhausted
            else
                cur[i] = 0;         // otherwise, reset this to 0 and increase the next dimension
        }
        else
            return true;            // current dimension not yet exhausted, return current vector
    }
    GISMO_ERROR("Something went wrong in nextLexicographic, wrong input?");
}


template<class Vec>
bool nextLexicographic(Vec& cur, const Vec& start, const Vec& end)
{
    const index_t d = cur.size();
    GISMO_ASSERT( d == start.size() && d == end.size(),
                  "Vector sizes don't match in nextLexicographic");

    for (index_t i = 0; i < d; ++i)
    {
        // increase current dimension
        if (++cur[i] == end[i])     // current dimension exhausted ?
        {
            if (i == d - 1)         // was it the last one?
                return false;       // then all elements exhausted
            else
                cur[i] = start[i];  // otherwise, reset this and increase the next dimension
        }
        else
            return true;            // current dimension not yet exhausted, return current vector
    }
    GISMO_ERROR("Something went wrong in nextLexicographic, wrong input?");
}


template<class Vec>
bool nextCubeVertex(Vec& cur, const Vec& start, const Vec& end)
{
    const int d = cur.size();
    GISMO_ASSERT( d == start.size() && d == end.size(),
                  "Vector sizes don't match in nextCubeVertex");

    for (int i = 0; i != d; ++i)
    {
        if ( cur[i] != end[i] )
        {
            cur[i] = end[i];
            return true;
        }
        else
            cur[i] = start[i];
    }
    return false;
}

template<class Vec>
bool nextCubeVertex(Vec& cur, const Vec& end)
{
    const index_t d = cur.size();
    GISMO_ASSERT( d == end.size(),
                  "Vector sizes don't match in nextCubeVertex");

    for (index_t i = 0; i != d; ++i)
    {
        if ( cur[i] != end[i] )
        {
            cur[i] = end[i];
            return true;
        }
        else
            cur[i] = 0;
    }
    return false;
}

template<class Vec>
bool nextCubeVertex(Vec& cur)
{
    const index_t d = cur.size();
    GISMO_ASSERT( (cur.array() >= 0).all() && (cur.array() <= 1).all(),
                  "Input must be a vector of zeros and ones, got: "<<cur.transpose() );

    for (index_t i = 0; i != d; ++i)
    {
        if ( cur[i] == 0 )
        {
            ++cur[i];
            return true;
        }
        else
            cur[i] = 0;
    }
    return false;
}

template<class Vec>
bool nextCubePoint(Vec& cur, const Vec& end)
{
    const index_t d = cur.size();
    GISMO_ASSERT(d == static_cast<int>(end.size()),
                 "Vector sizes don't match in nextCubePoint");

    for (index_t i = 0; i != d; ++i)
    {
        if ( cur[i] != end[i] )
        {
            ++cur[i];
            return true;
        }
        else
            cur[i] = 0;
    }
    --cur[0];//stopping flag
    return false;
}

template<class Vec>
bool nextCubePoint(Vec& cur, const Vec& start, const Vec& end)
{
    const index_t d = cur.size();
    GISMO_ASSERT( d == static_cast<int>(start.size()) &&
                  d == static_cast<int>(end.size()),
                  "Vector sizes don't match in nextCubePoint");

    for (index_t i = 0; i != d; ++i)
    {
        if ( cur[i] != end[i] )
        {
            ++cur[i];
            return true;
        }
        else
            cur[i] = start[i];

    }
    --cur[0];//stopping flag
    return false;
}

template<class Vec>
bool nextCubePoint(Vec & ci, const Vec& start, const Vec& end, const Vec& str)
{
    const index_t d = end.size();
    GISMO_ASSERT( d == static_cast<int>(start.size()) &&
                  d == static_cast<int>(end.size()),
                  "Vector sizes don't match in nextCubePoint");

    for (index_t i = 0; i != d; ++i)
    {
        if ( ci % str[i] != end[i] ) // cur[i] != end[i]
        {
            ci += str[i]; //++cur[i]
            return true;
        }
        else
            ci -= (end[i]-start[i])*str[i]; // cur[i] = start[i]
    }
    ci = end.dot(str)+1;//stopping flag
    return false;
}

template<class Vec>
bool nextCubeBoundary(Vec& cur, const Vec& start, const Vec& end)
{
    const index_t d = cur.size();
    GISMO_ASSERT( d == start.size() && d == end.size(),
                  "Vector sizes don't match in nextCubeBoundary");

    for (index_t i = 0; i != d; ++i)
    {        
        if ( cur[i] != end[i] )
        {
            if ( cur[i] == start[i] && ( i!=d-1 || d==1) ) // || d==1 to treat 1D
            {
                int c=i+1;
                for (int j = c; j!=d; ++j)
                    if ( (cur[j] == start[j]) || 
                         (cur[j] == end[j]) )
                        c++;
                
                if ( c==1 )
                    cur[i] = end[i];
                else
                    cur[i]++;
            }
            else
                cur[i]++;

            for (int k = i-1; k!=-1; --k)
                cur[k] = start[k];
            return true;
        }
    }
    return false;
}

template<class Vec>
bool nextCubeBoundaryOffset(Vec& cur, const Vec& start, const Vec& end, Vec & offset)
{
    const index_t d = cur.size();
    GISMO_ASSERT( d == start.size() && d == end.size(),
                  "Vector sizes don't match in nextCubeBoundaryOffset");

    for (index_t i = 0; i != d; ++i)
    {        
        if ( cur[i] != end[i] )
        {
            if ( cur[i] == start[i]+offset[i] && ( i!=d-1 || d==1) ) // || d==1 to treat 1D
            {
                int c = i+1;
                for (int j = c; j!=d; ++j)
                    if (cur[j] <=  start[j] + offset[j] ||
                        cur[j] >=    end[j] - offset[j] )
                        c++;
                
                if ( c==1 )
                    cur[i] = end[i] - offset[i];
                else
                    cur[i]++;
            }
            else
                cur[i]++;

            for (int k = i-1; k!=-1; --k)
                cur[k] = start[k];
            return true;
        }
    }

    return false;
}

template<class Vec>
bool nextCubeBoundaryOffset(Vec& cur, const Vec& start, const Vec& end, 
                            Vec & loffset, Vec & uoffset)
{
    const index_t d = cur.size();
    GISMO_ASSERT( d == start.size() && d == end.size(),
                  "Vector sizes don't match in nextCubeBoundaryOffset");

    for (index_t i = 0; i != d; ++i)
    {        
        if ( cur[i] != end[i] )
        {
            if ( cur[i] == start[i]+loffset[i] && ( i!=d-1 || d==1) ) // || d==1 to treat 1D
            {
                int c = i+1;
                for (int j = c; j!=d; ++j)
                    if (cur[j] <=  start[j] + loffset[j] ||
                        cur[j] >=    end[j] - uoffset[j] )
                        c++;
                
                if ( c==1 )
                    cur[i] = end[i] - uoffset[i];
                else
                    cur[i]++;
            }
            else
                cur[i]++;

            for (int k = i-1; k!=-1; --k)
                cur[k] = start[k];
            return true;
        }
    }

    return false;
}

inline index_t numCubeElements(const index_t k, const index_t d)
{
    GISMO_ASSERT(k >= 0 && k<=d, "Invalid arguments.");
    return binomial(d,k) * (1<<(d-k));
}

template<class Vec>
inline index_t dimCubeElement(const Vec & cur)
{
    return (cur.array() == 2).count();
}

template<class Vec>
void firstCubeElement(Vec & cur, const index_t k = 0)
{
    const index_t d = cur.size();
    GISMO_ASSERT(k >= 0 && k<=d, "Invalid arguments.");
    cur.topRows   (k  ).setConstant(2);
    cur.bottomRows(d-k).setConstant(0);
}

template<class Vec>
bool nextCubeElement(Vec & cur, const index_t k)
{
    const index_t d = cur.size();
    GISMO_ASSERT(k >= 0 && k<=cur.size(), "Invalid arguments.");

    index_t i;
    do
    {
        for (i = 0; i != d; ++i)
        {
            if ( cur[i] != 2 )
            {
                ++cur[i];
                if ( (cur.array() == 2).count() == k ) // dimCubeElement(cur)==k ?
                    return true;
                else
                    break;// skip face 
            }
            else
                cur[i] = 0;
        }
    }
    while ( i!=d );

    return false;
}

template <typename Z, int d>
void cubeIsometry( const gsVector<bool,d>    & flip,
                   const gsVector<index_t,d> & perm, 
                   gsVector<Z> & result)
{
    const index_t dd = flip.size(); //binary sequence of length d
    GISMO_ASSERT( dd == perm.size(), "Dimensions do not match in cubeIsometry");
    GISMO_ASSERT( perm.sum() == dd*(dd-1)/2, "Error in the permutation: "<< perm.transpose());

    gsVector<index_t,d> pstr(dd);
    for (index_t k=0; k!=dd; ++k)
        pstr[k] = (1<<perm[k]);

    result.resize(1<<dd);
    index_t r  = 0;
    index_t i;
    gsVector<bool,d> v = gsVector<bool,d>::Zero(dd);
    do
    {
        Z & c = result[r++];
        c = 0;
        for (index_t k=0; k!=dd; ++k)
            c += ( flip[perm[k]] == v[k] ) * pstr[k];
        
        for (i = 0; i != dd; ++i)
        {
            if ( !v[i] )
            {
                v[i] = true;
                break;//for
            }
            else
                v[i] = false;
        }
    }
    while (i!=dd);
}

template <int d>
void cubeIsometryMatrix ( const gsVector<bool,d>    & flip,
                          const gsVector<index_t,d> & perm,
                          gsMatrix<int,d,d> & result)
{
    result.setZero(d,d);
    for(index_t i = 0; i < d; ++i)
        result(perm(i),i) = (flip(perm(i)) ? 1 : -1);
}

template<class Vec>
void firstComposition( typename Vec::Scalar sum, index_t dim, Vec & res)
{  
    res.derived().resize(dim);
    res.setZero();
    res[0] = sum;
}

template<class Vec>
inline bool nextComposition(Vec & v)
{
    const index_t k = v.size() - 1;
    
    if (v[k] != v.sum())
    {
        for (index_t i = 0; i <= k; ++i)
        {
            if ( v[i]!=0 )
            {
                const typename Vec::Scalar t = v[i];
                v[i]    = 0  ;
                v[0]    = t-1;
                v[i+1] += 1  ;
                return true;
            }
        }
    }
    return false;
}


inline unsigned numCompositions(int sum, short_t dim)
{
    return binomial(sum+dim-1,dim-1);
}


template<class Vec, class Mat>
void firstMultiComposition(const Vec & a, index_t k, Mat & res)
{  
    const index_t d = a.size();
    res.setZero(k, d);
    res.row(0) = a.transpose();
}

template<class Mat>
inline bool nextMultiComposition(Mat & m)
{
    const index_t k = m.rows();
    const index_t d = m.cols();
           
    for (index_t j = 0; j != d; ++j)
    {
        typename Mat::ColXpr c_j = m.col(j);
        if ( nextComposition(c_j) )
            return true;
        else
        {
            const index_t n_j = c_j.sum();
            firstComposition(n_j, k, c_j);
        }
    }
    return false;
}


template<class Vec>
unsigned numMultiCompositions(const Vec & a, index_t k)
{
    unsigned result = 1;
    const index_t d = a.size();
    for (index_t j = 0; j != d; ++j)
        result *= binomial<unsigned>(a[j]+k-1, k-1);
    return result;
}


// Lexicographic comparison for gsVectors
template<class Z, unsigned d>
struct lex_less
{
    bool operator() (const gsVector<Z,d> & lhs, const gsVector<Z,d> & rhs) const
    {
        unsigned i = 0;
        while( (i<d) && (lhs[i] == rhs[i++]) ) ;
        return ( --i==d ? false : lhs[i] < rhs[i] );
    }
};


} // namespace gismo