Combination of a Collection with Repetitions

Question

There are a lot of links on http://stackoverflow.com for how to do combinations: Generating combinations in c++ But these links presume to draw from an infinite set without repetition. When given a finite collection which does have repetition, these algorithms construct duplicates. For example you can see the accepted solution to the linked question failing on a test case I constructed here: http://ideone.com/M7CyFc

Given the input set: vector<int> row {40, 40, 40, 50, 50, 60, 100};

I expect to see:

• 40 40 40
• 40 40 50
• 40 40 60
• 40 40 100
• 40 50 50
• 40 50 60
• 40 50 100
• 40 60 100
• 50 50 60
• 50 50 100
• 50 60 100

Obviously I can use the old method store the output and check for duplicates as I generate, but this requires a lot of extra memory and processing power. Is there an alternative that C++ provides me?

Answer 1

Combinations by definition do not respect order. This frees us to arrange the numbers in any order we see fit. Most notably we can rely on to provide a combination rank. Certainly the most logical way to rank combinations is in sorted order, so we'll be depending upon our inputs being sorted.

There is already precedent for this in the standard library. For example lower_bound which we will actually use in this solution. When used generally this may however require the user to sort before passing.

The function we will write to do this will take in iterators to the sorted collection which the next combination is to be drawn from, and iterators to the current combination. We'd also need the size but that can be derived from the distance between the combination's iterators.

template <typename InputIt, typename OutputIt>
bool next_combination(InputIt inFirst, InputIt inLast, OutputIt outFirst, OutputIt outLast) {
    assert(distance(inFirst, inLast) >= distance(outFirst, outLast));

    const auto front = make_reverse_iterator(outFirst);
    const auto back = make_reverse_iterator(outLast);
    auto it = mismatch(back, front, make_reverse_iterator(inLast)).first;

    const auto result = it != front;

    if (result) {
        auto ub = upper_bound(inFirst, inLast, *it);

        copy(ub, next(ub, distance(back, it) + 1), next(it).base());
    }
    return result;
}

This function is written in the format of the other algorithm functions, so any container that supports bidirectional iterators can be used with it. For our example though we'll use: const vector<unsigned int> row{ 40U, 40U, 40U, 50U, 50U, 60U, 100U }; which is, necessarily, sorted:

vector<unsigned int> it{ row.cbegin(), next(row.cbegin(), 3) };

do {
    copy(it.cbegin(), it.cend(), ostream_iterator<unsigned int>(cout, " "));
    cout << endl;
} while(next_combination(row.cbegin(), row.cend(), it.begin(), it.end()));

Live Example

---

After writing this answer I've done a bit more research and found N2639 which proposes a standardized next_combination, which was:

• Actively under consideration for a future TR, when work on TR2 was deferred pending
• Viewed positively at the time
• Due at least one more revision before any adoption
• Needed some reworking to reflect the addition of C++11 language facilities

[Source]

Using N2639's reference implementation requires mutability, so we'll use: vector<unsigned int> row{ 40U, 40U, 40U, 50U, 50U, 60U, 100U };. And our example code becomes:

vector<unsigned int>::iterator it = next(row.begin(), 3);

do {
    copy(row.begin(), it, ostream_iterator<unsigned int>(cout, " "));
    cout << endl;
} while(next_combination(row.begin(), it, row.end()));

Live Example

Answer 2

You can do something like this (maybe avoiding the recursion):

#include <iostream>
#include <vector>
#include <algorithm>

using std::cout;
using std::vector;

void perm( const vector<int> &v, vector<vector<int>> &p, vector<int> &t, int k, int d) {
    for ( int i = k; i < v.size(); ++i ) {
        // skip the repeted value
        if ( i != k && v[i] == v[i-1]) continue;
        t.push_back(v[i]);
        if ( d > 0 ) perm(v,p,t,i+1,d-1);
        else p.push_back(t);
        t.pop_back();
    }
}

int main() {
    int r = 3;
    vector<int> row {40, 40, 40, 50, 50, 60, 100};
    vector<vector<int>> pp;
    vector<int> pe;

    std::sort(row.begin(),row.end()); // that's necessary
    perm(row,pp,pe,0,r-1);

    cout << pp.size() << '\n';
    for ( auto & v : pp ) {
        for ( int i : v ) {
            cout << ' ' << i;
        }
        cout << '\n';
    }
    return 0;
}

Which outputs:

11
 40 40 40
 40 40 50
 40 40 60
 40 40 100
 40 50 50
 40 50 60
 40 50 100
 40 60 100
 50 50 60
 50 50 100
 50 60 100

I know, it's far from efficient, but if you get the idea you may come out with a better implementation.

Answer 3

Here is a class I once wrote in my university times to handle bosons. It's quite long, but it's generally usable and seems to work well. Additionally, it also gives ranking and unranking functionality. Hope that helps -- but don't ever ask me what I was doing back then ... ;-)

struct SymmetricIndex
{
    using StateType = std::vector<int>;
    using IntegerType = int;
    int M;
    int N;
    StateType Nmax;
    StateType Nmin;
    IntegerType _size;
    std::vector<IntegerType> store;
    StateType state;
    IntegerType _index;

    SymmetricIndex() = default;
    SymmetricIndex(int _M, int _N, int _Nmax = std::numeric_limits<int>::max(), int _Nmin = 0)
        : SymmetricIndex(_M, _N, std::vector<int>(_M + 1, std::min(_Nmax, _N)), StateType(_M + 1, std::max(_Nmin, 0)))
    {}

    SymmetricIndex(int _M, int _N, StateType const& _Nmax, StateType const& _Nmin)
        : N(_N)
        , M(_M)
        , Nmax(_Nmax)
        , Nmin(_Nmin)
        , store(addressArray())
        , state(M)
        , _index(0)
    {
        reset();
        _size = W(M, N);
    }

    friend std::ostream& operator<<(std::ostream& os, SymmetricIndex const& sym);

    SymmetricIndex& reset()
    {
        return setBegin();
    }


    bool setBegin(StateType& state, StateType const& Nmax, StateType const& Nmin) const
    {
        int n = N;
        for (int i = 0; i<M; ++i)
        {
            state[i] = Nmin[i];
            n -= Nmin[i];
        }

        for (int i = 0; i<M; ++i)
        {
            state[i] = std::min(n + Nmin[i], Nmax[i]);
            n -= Nmax[i] - Nmin[i];
            if (n <= 0)
                break;
        }

        return true;
    }


    SymmetricIndex& setBegin()
    {
        setBegin(state, Nmax, Nmin);
        _index = 0;
        return *this;
    }


    bool isBegin() const
    {
        return _index==0;
    }

    bool setEnd(StateType& state, StateType const& Nmax, StateType const& Nmin) const
    {
        int n = N;
        for (int i = 0; i < M; ++i)
        {
            state[i] = Nmin[i];
            n -= Nmin[i];
        }

        for (int i = M - 1; i >= 0; --i)
        {
            state[i] = std::min(n + Nmin[i], Nmax[i]);
            n -= Nmax[i] - Nmin[i];
            if (n <= 0)
                break;
        }
        return true;
    }
    SymmetricIndex& setEnd()
    {
        setEnd(state, Nmax, Nmin);
        _index = _size - 1;
        return *this;
    }

    bool isEnd() const
    {
        return _index == _size-1;
    }

    IntegerType index() const
    {
        return _index;
    }

    IntegerType rank(StateType const& state) const
    {
        IntegerType ret = 0;
        int n = 0;

        for (int i = 0; i < M; ++i)
        {
            n += state[i];
            for (int k = Nmin[i]; k < state[i]; ++k)
                ret += store[(n - k) * M + i];
        }

        return ret;
    }

    IntegerType rank() const
    {
        return rank(state);
    }

    StateType unrank(IntegerType rank) const
    {
        StateType ret(M);

        int n = N;
        for (int i = M-1; i >= 0; --i)
        {
            int ad = 0;
            int k = std::min(Nmax[i] - 1, n);

            for (int j = Nmin[i]; j <= k; ++j)
                ad+=store[(n - j) * M + i];

            while (ad > rank && k >= Nmin[i])
            {
                ad -= store[(n - k) * M + i];
                --k;
            }

            rank -= ad;
            ret[i] = k+1;
            n -= ret[i];
            if (n <= 0)
            {
                return ret;
            }
        }

        return ret;
    }

    IntegerType size() const
    {
        return _size;
    }

    operator StateType& ()
    {
        return state;
    }

    auto operator[](int i) -> StateType::value_type& { return state[i]; }

    operator StateType const& () const
    {
        return state;
    }
    auto operator[](int i) const -> StateType::value_type const& { return state[i]; }


    bool nextState(StateType& state, StateType const& Nmax, StateType const& Nmin) const
    {
        //co-lexicographical ordering with Nmin and Nmax:

        // (1) find first position which can be decreased
        //     then we have state[k] = Nmin[k]  for k in [0,pos]
        int pos = M - 1;
        for (int k = 0; k < M - 1; ++k)
        {
            if (state[k] > Nmin[k])
            {
                pos = k;
                break;
            }
        }

        // if nothing found to decrease, return
        if (pos == M - 1)
        {
            return false;
        }

        // (2) find first position after pos which can be increased
        //     then we have state[k] = Nmin[k]  for k in [0,pos]
        int next = 0;
        for (int k = pos + 1; k < M; ++k)
        {
            if (state[k] < Nmax[k])
            {
                next = k;
                break;
            }
        }
        if (next == 0)
        {
            return false;
        }
        --state[pos];
        ++state[next];

        // (3) get occupation in [pos,next-1] and set to Nmin[k]
        int n = 0;
        for (int k = pos; k < next; ++k)
        {
            n += state[k] - Nmin[k];
            state[k] = Nmin[k];
        }

        // (4) fill up from the start
        for (int i = 0; i<M; ++i)
        {
            if (n <= 0)
                break;
            int add = std::min(n, Nmax[i] - state[i]);
            state[i] += add;
            n -= add;
        }

        return true;
    }


    SymmetricIndex& operator++()
    {
        bool inc = nextState(state, Nmax, Nmin);
        if (inc) ++_index;
        return *this;
    }
    SymmetricIndex operator++(int)
    {
        auto ret = *this;
        this->operator++();
        return ret;
    }

    bool previousState(StateType& state, StateType const& Nmax, StateType const& Nmin) const
    {
        ////co-lexicographical ordering with Nmin and Nmax:

        // (1) find first position which can be increased
        //     then we have state[k] = Nmax[k]  for k in [0,pos-1]
        int pos = M - 1;
        for (int k = 0; k < M - 1; ++k)
        {
            if (state[k] < Nmax[k])
            {
                pos = k;
                break;
            }
        }

        // if nothing found to increase, return
        if (pos == M - 1)
        {
            return false;
        }

        // (2) find first position after pos which can be decreased
        //     then we have state[k] = Nmin[k]  for k in [pos+1,next]
        int next = 0;
        for (int k = pos + 1; k < M; ++k)
        {
            if (state[k] > Nmin[k])
            {
                next = k;
                break;
            }
        }
        if (next == 0)
        {
            return false;
        }
        ++state[pos];
        --state[next];

        int n = 0;
        for (int k = 0; k <= pos; ++k)
        {
            n += state[k] - Nmin[k];
            state[k] = Nmin[k];
        }
        if (n == 0)
        {
            return true;
        }

        for (int i = next-1; i>=0; --i)
        {
            int add = std::min(n, Nmax[i] - state[i]);
            state[i] += add;
            n -= add;
            if (n <= 0)
                break;
        }

        return true;
    }


    SymmetricIndex operator--()
    {
        bool dec = previousState(state, Nmax, Nmin);
        if (dec) --_index;
        return *this;
    }
    SymmetricIndex operator--(int)
    {
        auto ret = *this;
        this->operator--();
        return ret;
    }

    int multinomial() const
    {
        auto v = const_cast<std::remove_reference<decltype(state)>::type&>(state);
        return multinomial(v);
    }

    int multinomial(StateType& state) const
    {
        int ret = 1;
        int n = state[0];
        for (int i = 1; i < M; ++i)
        {
            n += state[i];
            ret *= binomial(n, state[i]);
        }
        return ret;
    }

    SymmetricIndex& random(StateType const& _Nmin)
    {
        static std::mt19937 rng;

        state = _Nmin;
        int n = std::accumulate(std::begin(state), std::end(state), 0);
        auto weight = [&](int i) { return state[i] < Nmax[i] ? 1 : 0; };

        for (int i = n; i < N; ++i)
        {
            std::discrete_distribution<int> d(N, 0, N, weight);
            ++state[d(rng)];
        }
        _index = rank();

        return *this;
    }
    SymmetricIndex& random()
    {
        return random(Nmin);
    }

private:
    IntegerType W(int m, int n) const
    {
        if (m < 0 || n < 0) return 0;
        else if (m == 0 && n == 0) return 1;
        else if (m == 0 && n > 0) return 0;
        //else if (m > 0 && n < Nmin[m-1]) return 0;
        else
        {
            //static std::map<std::tuple<int, int>, IntegerType> memo;

            //auto it = memo.find(std::make_tuple(k, m));
            //if (it != std::end(memo))
            //{
            //  return it->second;
            //}

            IntegerType ret = 0;
            for (int i = Nmin[m-1]; i <= std::min(Nmax[m-1], n); ++i)
                ret += W(m - 1, n - i);

            //memo[std::make_tuple(k, m)] = ret;
            return ret;
        }
    }

    IntegerType binomial(int m, int n) const
    {
        static std::vector<int> store;
        if (store.empty())
        {
            std::function<IntegerType(int, int)> bin = [](int n, int k)
            {
                int res = 1;

                if (k > n - k)
                    k = n - k;

                for (int i = 0; i < k; ++i)
                {
                    res *= (n - i);
                    res /= (i + 1);
                }

                return res;
            };

            store.resize(M*M);
            for (int i = 0; i < M; ++i)
            {
                for (int j = 0; j < M; ++j)
                {
                    store[i*M + j] = bin(i, j);
                }
            }
        }
        return store[m*M + n];
    }

    auto addressArray() const -> std::vector<int>
    {
        std::vector<int> ret((N + 1) * M);

        for (int n = 0; n <= N; ++n)
        {
            for (int m = 0; m < M; ++m)
            {
                ret[n*M + m] = W(m, n);
            }
        }
        return ret;
    }
};

std::ostream& operator<<(std::ostream& os, SymmetricIndex const& sym)
{
    for (auto const& i : sym.state)
    {
        os << i << " ";
    }
    return os;
}

Use it like

int main()
{
    int M=4;
    int N=3;
    std::vector<int> Nmax(M, N);
    std::vector<int> Nmin(M, 0);
    Nmax[0]=3;
    Nmax[1]=2;
    Nmax[2]=1;
    Nmax[3]=1;

    SymmetricIndex sym(M, N, Nmax, Nmin);

    while(!sym.isEnd())
    {
        std::cout<<sym<<"   "<<sym.rank()<<std::endl;
        ++sym;
    }
    std::cout<<sym<<"   "<<sym.rank()<<std::endl;
}

This will output

3 0 0 0    0    (corresponds to {40,40,40})
2 1 0 0    1    (-> {40,40,50})
1 2 0 0    2    (-> {40,50,50})
2 0 1 0    3    ...
1 1 1 0    4
0 2 1 0    5
2 0 0 1    6
1 1 0 1    7
0 2 0 1    8
1 0 1 1    9
0 1 1 1    10   (-> {50,60,100})

DEMO

Note that I assumed here an ascending mapping of your set elements (i.e. the number 40's is given by index 0, the number of 50's by index 1, and so on).

---

---

More precisely: Turn your list into a map<std::vector<int>, int> like

std::vector<int> v{40,40,40,50,50,60,100};

std::map<int, int> m;

for(auto i : v)
{
    ++m[i];
}

Then use

int N = 3;
int M = m.size();
std::vector<int> Nmin(M,0);
std::vector<int> Nmax;
std::vector<int> val;

for(auto i : m)
{
    Nmax.push_back(m.second);
    val.push_back(m.first);
}

SymmetricIndex sym(M, N, Nmax, Nmin);

as input to the SymmetricIndex class.

To print the output, use

    while(!sym.isEnd())
    {
         for(int i=0; i<M; ++i)
         {
              for(int j = 0; j<sym[i]; ++j)
              {
                  std::cout<<val[i]<<" ";
              }
         }
         std::cout<<std::endl;
    }

    for(int i=0; i<M; ++i)
    {
        for(int j = 0; j<sym[i]; ++j)
        {
            std::cout<<val[i]<<" ";
        }
    }
    std::cout<<std::endl;

All untested, but it should give the idea.