IntegerComposition.hpp

Go to the documentation of this file.

#pragma once
#include "Headers/mydebug.hpp"
#include <Eigen/Dense>
#include <iostream>
 
#ifndef __NVCC__
    #define __host__
    #define __device__
#endif
 
#if __has_include(<omp.h>)
    #include <omp.h>
__host__ __device__ static inline int get_max_threads() {
    #ifdef __CUDA_ARCH__
    return 1;
    #else
    return omp_get_max_threads();
    #endif
}
__host__ __device__ static inline int get_thread_num() {
    #ifdef __CUDA_ARCH__
    return 0;
    #else
    return omp_get_thread_num();
    #endif
}
#else
constexpr static inline int get_max_threads() { return 1; }
constexpr static inline int get_thread_num() { return 0; }
#endif
 
template<class T>
class is_container {
        template<class U, void* dummy = (&U::operator[], &U::resize, nullptr)>
        static std::true_type  test(U*);
        static std::false_type test(...);
 
    public:
        static constexpr bool value = decltype(test((T*)nullptr))::value;
};
template<class T>
inline constexpr bool is_container_v = is_container<T>::value;
 
class IntegerComposition {
    private:
        typedef unsigned long   Index_t;
        int                     m_N;
        int                     m_Length;
        int                     m_Max;
        Index_t                 m_dim;
        Eigen::MatrixX<Index_t> workA;
        Eigen::MatrixX<Index_t> workB;
 
    public:
        __host__ __device__ IntegerComposition() { debug_constructor_printf(Default); }
        __host__ __device__ IntegerComposition(int N, int Length, int Max);
        __host__ __device__ IntegerComposition(IntegerComposition const& other)
            : m_N(other.value()),
              m_Length(other.length()),
              m_Max(other.max()),
              m_dim(other.dim()),
              workA(other.workA),
              workB(other.workB) {
            debug_constructor_printf((Copy));
        }
 
        __host__ __device__ int     value() const { return m_N; };
        __host__ __device__ int     length() const { return m_Length; };
        __host__ __device__ int     max() const { return m_Max; };
        __host__ __device__ Index_t dim() const { return m_dim; };
 
        __host__ void status(void) const {
            std::cout << "IntegerComposition.status()"
                      << "\n"
                      << "\tValue:" << this->m_N << "\n"
                      << "\tLength:" << this->m_Length << "\n"
                      << "\tMax:   " << this->m_Max << "\n"
                      << "\tDim:   " << this->m_dim << std::endl;
            std::cout << "workA:\n"
                      << workA << "\n"
                      << "workB:\n"
                      << workB << std::endl;
        }
 
        template<typename Vector_t>
        __host__ __device__ Index_t partToOrdinal(Vector_t const& vec) const;
 
        template<class Vector_t>
        __host__ __device__ void ordinalToPart(Vector_t& vec, Index_t const num) const;
 
        __host__ __device__ Eigen::RowVectorXi ordinalToPart(Index_t const num) const {
            Eigen::RowVectorXi res(this->m_Length);
            this->ordinalToPart(res, num);
            return res;
        }
 
        template<typename Index_>
        __host__ __device__ int locNumber(Index_ const num, Index_ const pos) const;
 
        template<class Vector_t>
        __host__ __device__ Index_t translate(Index_t const num, int trans,
                                              Vector_t& workSpace) const;
 
        __host__ __device__ Index_t translate(Index_t const num, int trans) const {
            Eigen::ArrayXi config(m_Length);
            return this->translate(num, trans, config);
        }
};
 
inline IntegerComposition::IntegerComposition(int N, int Length, int Max)
    : m_N(N),
      m_Length(Length),
      m_Max(Max < N ? Max : N),
      workA(N + 1, Length),
      workB(N + 1, Length) {
    // debug_printf("IntegerComposition::IntegerComposition(%d,%d,%d) : Dim = %d\n", (int)N,
    //              (int)Length, (int)Max, (int)m_dim);
    if(Length <= 1) {
        m_dim = Length;
        return;
    }
 
    int   max  = std::min(Max, N);
    auto& Dims = workB;
    for(int l = 0; l < Length; ++l) Dims(0, l) = 1;
    for(int n = N; n > max; --n) {
        Dims(n, 1)  = 0;
        Dims(n, 0)  = 0;
        workA(n, 0) = 0;
    }
    for(int n = max; n >= 0; --n) {
        Dims(n, 1)  = 1;
        Dims(n, 0)  = 0;
        workA(n, 0) = 0;
    }
    // Recursively calculate Dims(l,n) for remaining (l,n)
    for(int l = 2; l < Length; ++l)
        for(int n = 1; n <= N; ++n) {
            Dims(n, l) = 0;
            for(int k = 0; k <= std::min(max, n); ++k) { Dims(n, l) += Dims(n - k, l - 1); }
        }
    //   std::cout << Dims << "\n" << std::endl;
 
    for(int l = 1; l < Length; ++l) {
        workA(0, l) = Dims(N, Length - l);
        for(int n = 1; n <= N; ++n) workA(n, l) = workA(n - 1, l) + Dims(N - n, Length - l);
    }
    // std::cout << workA << "\n" << std::endl;
 
    for(int l = 0; l < Length; ++l) workB(0, l) = 0;
    for(int n = 1; n <= N; ++n) {
        for(int l = 1; l < Length - 1; ++l) workB(n, l) = workA(n - 1, l) - workA(n - 1, l + 1);
        workB(n, Length - 1) = workA(n - 1, Length - 1);
    }
    // std::cout << workB << "\n" << std::endl;
 
    m_dim = workA(max, 1);
    // std::cout << "Dim = " << m_dim << std::endl;
    debug_printf("(END) IntegerComposition::IntegerComposition(%d,%d,%d) : Dim = %d\n", (int)N,
                 (int)Length, (int)Max, (int)m_dim);
}
 
template<typename Vector_t>
__host__ __device__ inline IntegerComposition::Index_t IntegerComposition::partToOrdinal(
    Vector_t const& vec) const {
    IntegerComposition::Index_t z = 0, res = 0;
    for(int l = 1; l != int(m_Length); ++l) {
        assert(m_Length - l >= 0);
        assert(m_Length - l < vec.size());
        z += vec[m_Length - l];
 
        assert(z < workB.rows());
        assert(l < workB.cols());
        res += workB(z, l);
    }
    return res;
}
 
template<typename Vector_t>
__host__ __device__ inline void IntegerComposition::ordinalToPart(
    Vector_t& vec, IntegerComposition::Index_t const num) const {
    IntegerComposition::Index_t num_copy = num;
    int                         z = 0, zprev = 0;
#ifndef __CUDA_ARCH__
    vec.resize(m_Length);
#else
    assert(vec.size() >= m_Length);
#endif
    // if constexpr(is_container_v<Vector_t>) { vec.resize(m_Length); }
    for(int l = 1; l != m_Length; ++l) {
        assert(z < workA.rows());
        assert(l < workA.cols());
        while(workA(z, l) <= num_copy) z += 1;
 
        assert(m_Length - l >= 0);
        assert(m_Length - l < vec.size());
        vec[m_Length - l] = z - zprev;
        zprev             = z;
 
        assert(z < workB.rows());
        assert(l < workB.cols());
        num_copy -= workB(z, l);
    }
    vec[0] = m_N - z;
}
 
template<typename Index_>
__host__ __device__ inline int IntegerComposition::locNumber(Index_ const num,
                                                             Index_ const pos) const {
    // Eigen::ArrayXi config(m_Length);
    // this->ordinalToPart(config, num);
    // return config(pos);
 
    auto num_copy = num;
    auto position = ((pos % m_Length) + m_Length) % m_Length;
    int  z = 0, zprev = 0;
    int  l;
    for(l = 1; l != m_Length - position; ++l) {
        assert(z < workA.rows());
        assert(l < workA.cols());
        while(workA(z, l) <= num_copy) z += 1;
        zprev = z;
 
        assert(z < workB.rows());
        assert(l < workB.cols());
        num_copy -= workB(z, l);
    }
    assert(z < workA.rows());
    assert(l < workA.cols());
    if(pos == 0) { return m_N - z; }
    else {
        while(workA(z, l) <= num_copy) z += 1;
        return z - zprev;
    }
}
 
template<class Vector_t>
__host__ __device__ inline IntegerComposition::Index_t IntegerComposition::translate(
    IntegerComposition::Index_t const num, int trans, Vector_t& workSpace) const {
    if(workSpace.size() < m_Length)
        debug_printf("%s\n\tworkSpace.size()=%d\n", __PRETTY_FUNCTION__, (int)workSpace.size());
    assert(workSpace.size() >= m_Length);
    // debug_printf("%s\n\t workSpace=%p, workSpace.size()=%d, m_Length=%d\n", __PRETTY_FUNCTION__,
    //                workSpace.data(), (int)workSpace.size(), (int)m_Length);
    this->ordinalToPart(workSpace, num);
    Index_t z = 0, res = 0;
    for(int l = 0; l < int(m_Length); ++l) {
        // if(!(z < workB.rows()) || !(l < workB.cols())) {
        //  debug_printf("\t (z=%d, l=%d), workB.rows()=%d, workB.cols()=%d\n", (int)z, (int)l,
        //               (int)workB.rows(), (int)workB.cols());
        // }
        assert(z < workB.rows());
        assert(l < workB.cols());
        res += workB(z, l);
 
        assert(int(m_Length - 1 - l + trans) % workSpace.size() < workSpace.size());
        z += workSpace[int(m_Length - 1 - l + trans) % workSpace.size()];
    }
    return res;
}