multiseq_selection.h

Go to the documentation of this file.

/***************************************************************************
 * include/stxxl/bits/parallel/multiseq_selection.h
 *
 * Functions to find elements of a certain global rank in multiple sorted
 * sequences. Also serves for splitting such sequence sets.
 *
 * Extracted from MCSTL - http://algo2.iti.uni-karlsruhe.de/singler/mcstl/
 *
 * Part of the STXXL. See http://stxxl.sourceforge.net
 *
 * Copyright (C) 2007 Johannes Singler <[email protected]>
 * Copyright (C) 2014 Timo Bingmann <[email protected]>
 *
 * Distributed under the Boost Software License, Version 1.0.
 * (See accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 **************************************************************************/

#ifndef STXXL_PARALLEL_MULTISEQ_SELECTION_HEADER
#define STXXL_PARALLEL_MULTISEQ_SELECTION_HEADER

#include <stxxl/bits/namespace.h>
#include <stxxl/bits/parallel/types.h>
#include <stxxl/bits/parallel/compiletime_settings.h>

#include <vector>
#include <queue>
#include <cstdlib>
#include <algorithm>

STXXL_BEGIN_NAMESPACE

namespace parallel {

//! Compare a pair of types lexcigraphically, ascending.
template <typename T1, typename T2, typename Comparator>
class lexicographic
 : public std::binary_function<std::pair<T1, T2>, std::pair<T1, T2>, bool>
{
protected:
 Comparator& m_comp;

public:
 lexicographic(Comparator& comp) : m_comp(comp) { }

 inline bool operator () (const std::pair<T1, T2>& p1,
 const std::pair<T1, T2>& p2)
 {
 if (m_comp(p1.first, p2.first))
 return true;

 if (m_comp(p2.first, p1.first))
 return false;

 // firsts are equal
 return p1.second < p2.second;
 }
};

//! Compare a pair of types lexcigraphically, descending.
template <typename T1, typename T2, typename Comparator>
class lexicographic_rev
 : public std::binary_function<std::pair<T1, T2>, std::pair<T1, T2>, bool>
{
protected:
 Comparator& m_comp;

public:
 lexicographic_rev(Comparator& comp) : m_comp(comp) { }

 inline bool operator () (const std::pair<T1, T2>& p1,
 const std::pair<T1, T2>& p2)
 {
 if (m_comp(p2.first, p1.first))
 return true;

 if (m_comp(p1.first, p2.first))
 return false;

 // firsts are equal
 return p2.second < p1.second;
 }
};

/*!
 * Splits several sorted sequences at a certain global rank, resulting in a
 * splitting point for each sequence. The sequences are passed via a sequence
 * of random-access iterator pairs, none of the sequences may be empty. If
 * there are several equal elements across the split, the ones on the left side
 * will be chosen from sequences with smaller number.
 *
 * \param begin_seqs Begin of the sequence of iterator pairs.
 * \param end_seqs End of the sequence of iterator pairs.
 * \param rank The global rank to partition at.
 * \param begin_offsets A random-access sequence begin where the result will be
 * stored in. Each element of the sequence is an iterator that points to the
 * first element on the greater part of the respective sequence.
 * \param comp The ordering functor, defaults to std::less<T>.
 */
template <typename RanSeqs, typename RankType, typename RankIterator,
 typename Comparator>
void multiseq_partition(
 const RanSeqs& begin_seqs, const RanSeqs& end_seqs,
 const RankType& rank,
 RankIterator begin_offsets,
 Comparator comp = std::less<
 typename std::iterator_traits<typename std::iterator_traits<RanSeqs>
 ::value_type::first_type>::value_type
 >()) //std::less<T>
{
 STXXL_PARALLEL_PCALL(end_seqs - begin_seqs);

 typedef typename std::iterator_traits<RanSeqs>
 ::value_type::first_type iterator;
 typedef typename std::iterator_traits<iterator>
 ::difference_type diff_type;
 typedef typename std::iterator_traits<iterator>
 ::value_type value_type;

 typedef std::pair<value_type, diff_type> sample_pair;
 // comparators for sample_pair
 lexicographic<value_type, diff_type, Comparator> lcomp(comp);
 lexicographic_rev<value_type, diff_type, Comparator> lrcomp(comp);

 // number of sequences, number of elements in total (possibly including padding)
 const diff_type m = std::distance(begin_seqs, end_seqs);
 diff_type nmax, n;
 RankType N = 0;

 for (diff_type i = 0; i < m; ++i)
 {
 N += std::distance(begin_seqs[i].first, begin_seqs[i].second);
 assert(std::distance(begin_seqs[i].first, begin_seqs[i].second) > 0);
 }

 if (rank == N)
 {
 for (diff_type i = 0; i < m; ++i)
 begin_offsets[i] = begin_seqs[i].second; // very end
 return;
 }

 assert(m != 0 && N != 0 && rank < N);

 diff_type* seqlen = new diff_type[m];

 seqlen[0] = std::distance(begin_seqs[0].first, begin_seqs[0].second);
 nmax = seqlen[0];
 for (diff_type i = 1; i < m; ++i)
 {
 seqlen[i] = std::distance(begin_seqs[i].first, begin_seqs[i].second);
 nmax = std::max(nmax, seqlen[i]);
 }

 // pad all lists to this length, at least as long as any ns[i], equality
 // iff nmax = 2^k - 1
 diff_type l = round_up_to_power_of_two(nmax + 1) - 1;

 diff_type* a = new diff_type[m], * b = new diff_type[m];

 for (diff_type i = 0; i < m; ++i)
 a[i] = 0, b[i] = l;

 n = l / 2;

 // invariants:
 // 0 <= a[i] <= seqlen[i], 0 <= b[i] <= l

#define S(i) (begin_seqs[i].first)

 // initial partition

 std::vector<sample_pair> sample;

 for (diff_type i = 0; i < m; ++i) {
 if (n < seqlen[i]) {
 // sequence long enough
 sample.push_back(sample_pair(S(i)[n], i));
 }
 }

 std::sort(sample.begin(), sample.end(), lcomp);

 for (diff_type i = 0; i < m; ++i) {
 // conceptual infinity
 if (n >= seqlen[i]) {
 // sequence too short, conceptual infinity
 sample.push_back(sample_pair(S(i)[0] /*dummy element*/, i));
 }
 }

 diff_type localrank = rank / l;

 diff_type j;
 for (j = 0; j < localrank && ((n + 1) <= seqlen[sample[j].second]); ++j)
 a[sample[j].second] += n + 1;
 for ( ; j < m; ++j)
 b[sample[j].second] -= n + 1;

 // further refinement

 while (n > 0)
 {
 n /= 2;

 const value_type* lmax = NULL; // impossible to avoid the warning?
 for (diff_type i = 0; i < m; ++i)
 {
 if (a[i] > 0)
 {
 if (!lmax)
 lmax = &(S(i)[a[i] - 1]);
 else
 {
 // max, favor rear sequences
 if (!comp(S(i)[a[i] - 1], *lmax))
 lmax = &(S(i)[a[i] - 1]);
 }
 }
 }

 for (diff_type i = 0; i < m; ++i)
 {
 diff_type middle = (b[i] + a[i]) / 2;
 if (lmax && middle < seqlen[i] && comp(S(i)[middle], *lmax))
 a[i] = std::min(a[i] + n + 1, seqlen[i]);
 else
 b[i] -= n + 1;
 }

 diff_type leftsize = 0;
 for (diff_type i = 0; i < m; ++i)
 leftsize += a[i] / (n + 1);

 diff_type skew = rank / (n + 1) - leftsize;

 if (skew > 0)
 {
 // move to the left, find smallest
 std::priority_queue<sample_pair, std::vector<sample_pair>,
 lexicographic_rev<value_type, diff_type, Comparator> >
 pq(lrcomp);

 for (diff_type i = 0; i < m; ++i) {
 if (b[i] < seqlen[i])
 pq.push(sample_pair(S(i)[b[i]], i));
 }

 for ( ; skew != 0 && !pq.empty(); --skew)
 {
 diff_type source = pq.top().second;
 pq.pop();

 a[source] = std::min(a[source] + n + 1, seqlen[source]);
 b[source] += n + 1;

 if (b[source] < seqlen[source])
 pq.push(sample_pair(S(source)[b[source]], source));
 }
 }
 else if (skew < 0)
 {
 // move to the right, find greatest
 std::priority_queue<sample_pair, std::vector<sample_pair>,
 lexicographic<value_type, diff_type, Comparator> >
 pq(lcomp);

 for (diff_type i = 0; i < m; ++i) {
 if (a[i] > 0)
 pq.push(sample_pair(S(i)[a[i] - 1], i));
 }

 for ( ; skew != 0; ++skew)
 {
 diff_type source = pq.top().second;
 pq.pop();

 a[source] -= n + 1;
 b[source] -= n + 1;

 if (a[source] > 0)
 pq.push(sample_pair(S(source)[a[source] - 1], source));
 }
 }
 }

 // postconditions: a[i] == b[i] in most cases, except when a[i] has been
 // clamped because of having reached the boundary

 // now return the result, calculate the offset, compare the keys on both
 // edges of the border

 // maximum of left edge, minimum of right edge
 value_type* maxleft = NULL, * minright = NULL;
 for (diff_type i = 0; i < m; ++i)
 {
 if (a[i] > 0)
 {
 if (!maxleft)
 {
 maxleft = &(S(i)[a[i] - 1]);
 }
 else
 {
 // max, favor rear sequences
 if (!comp(S(i)[a[i] - 1], *maxleft))
 maxleft = &(S(i)[a[i] - 1]);
 }
 }
 if (b[i] < seqlen[i])
 {
 if (!minright)
 {
 minright = &(S(i)[b[i]]);
 }
 else
 {
 // min, favor fore sequences
 if (comp(S(i)[b[i]], *minright))
 minright = &(S(i)[b[i]]);
 }
 }
 }

 for (diff_type i = 0; i < m; ++i)
 begin_offsets[i] = S(i) + a[i];

 delete[] seqlen;
 delete[] a;
 delete[] b;
}

/*!
 * Selects the element at a certain global rank from several sorted sequences.
 *
 * The sequences are passed via a sequence of random-access iterator pairs,
 * none of the sequences may be empty.
 *
 * \param begin_seqs Begin of the sequence of iterator pairs.
 * \param end_seqs End of the sequence of iterator pairs.
 * \param rank The global rank to partition at.
 * \param offset The rank of the selected element in the global subsequence of
 * elements equal to the selected element. If the selected element is unique,
 * this number is 0.
 * \param comp The ordering functor, defaults to std::less. */
template <typename ValueType, typename RanSeqs, typename RankType, typename Comparator>
ValueType multiseq_selection(const RanSeqs& begin_seqs, const RanSeqs& end_seqs,
 const RankType& rank,
 RankType& offset,
 Comparator comp = std::less<ValueType>())
{
 STXXL_PARALLEL_PCALL(end_seqs - begin_seqs);

 typedef typename std::iterator_traits<RanSeqs>
 ::value_type::first_type iterator;
 typedef typename std::iterator_traits<iterator>
 ::difference_type diff_type;

 typedef std::pair<ValueType, diff_type> sample_pair;
 // comparators for sample_pair
 lexicographic<ValueType, diff_type, Comparator> lcomp(comp);
 lexicographic_rev<ValueType, diff_type, Comparator> lrcomp(comp);

 // number of sequences, number of elements in total (possibly including padding)
 const diff_type m = std::distance(begin_seqs, end_seqs);
 diff_type nmax, n;
 RankType N = 0;

 for (diff_type i = 0; i < m; ++i)
 N += std::distance(begin_seqs[i].first, begin_seqs[i].second);

 if (m == 0 || N == 0 || rank < 0 || rank >= N)
 // result undefined when there is no data or rank is outside bounds
 throw std::exception();

 diff_type* seqlen = new diff_type[m];

 seqlen[0] = std::distance(begin_seqs[0].first, begin_seqs[0].second);
 nmax = seqlen[0];
 for (diff_type i = 1; i < m; ++i)
 {
 seqlen[i] = std::distance(begin_seqs[i].first, begin_seqs[i].second);
 nmax = std::max(nmax, seqlen[i]);
 }

 // pad all lists to this length, at least as long as any ns[i], equliaty iff
 // nmax = 2^k - 1
 diff_type l = round_up_to_power_of_two(nmax + 1) - 1;

 diff_type* a = new diff_type[m], * b = new diff_type[m];

 for (diff_type i = 0; i < m; ++i)
 a[i] = 0, b[i] = l;

 n = l / 2;

 // invariants:
 // 0 <= a[i] <= seqlen[i], 0 <= b[i] <= l

#define S(i) (begin_seqs[i].first)

 //initial partition

 std::vector<sample_pair> sample;

 for (diff_type i = 0; i < m; ++i) {
 if (n < seqlen[i])
 sample.push_back(sample_pair(S(i)[n], i));
 }

 std::sort(sample.begin(), sample.end(), lcomp);

 for (diff_type i = 0; i < m; ++i) {
 //conceptual infinity
 if (n >= seqlen[i])
 sample.push_back(sample_pair(S(i)[0] /*dummy element*/, i));
 }

 diff_type localrank = rank / l;

 diff_type j;
 for (j = 0; j < localrank && ((n + 1) <= seqlen[sample[j].second]); ++j)
 a[sample[j].second] += n + 1;
 for ( ; j < m; ++j)
 b[sample[j].second] -= n + 1;

 // further refinement

 while (n > 0)
 {
 n /= 2;

 const ValueType* lmax = NULL;
 for (diff_type i = 0; i < m; ++i)
 {
 if (a[i] > 0)
 {
 if (!lmax)
 {
 lmax = &(S(i)[a[i] - 1]);
 }
 else
 {
 // max
 if (comp(*lmax, S(i)[a[i] - 1]))
 lmax = &(S(i)[a[i] - 1]);
 }
 }
 }

 for (diff_type i = 0; i < m; ++i)
 {
 diff_type middle = (b[i] + a[i]) / 2;
 if (lmax && middle < seqlen[i] && comp(S(i)[middle], *lmax))
 a[i] = std::min(a[i] + n + 1, seqlen[i]);
 else
 b[i] -= n + 1;
 }

 diff_type leftsize = 0;
 for (diff_type i = 0; i < m; ++i)
 leftsize += a[i] / (n + 1);

 diff_type skew = rank / (n + 1) - leftsize;

 if (skew > 0)
 {
 // move to the left, find smallest
 std::priority_queue<sample_pair, std::vector<sample_pair>,
 lexicographic_rev<ValueType, diff_type, Comparator> >
 pq(lrcomp);

 for (diff_type i = 0; i < m; ++i) {
 if (b[i] < seqlen[i])
 pq.push(sample_pair(S(i)[b[i]], i));
 }

 for ( ; skew != 0 && !pq.empty(); --skew)
 {
 diff_type source = pq.top().second;
 pq.pop();

 a[source] = std::min(a[source] + n + 1, seqlen[source]);
 b[source] += n + 1;

 if (b[source] < seqlen[source])
 pq.push(sample_pair(S(source)[b[source]], source));
 }
 }
 else if (skew < 0)
 {
 // move to the right, find greatest
 std::priority_queue<sample_pair, std::vector<sample_pair>,
 lexicographic<ValueType, diff_type, Comparator> >
 pq(lcomp);

 for (diff_type i = 0; i < m; ++i) {
 if (a[i] > 0)
 pq.push(sample_pair(S(i)[a[i] - 1], i));
 }

 for ( ; skew != 0; ++skew)
 {
 diff_type source = pq.top().second;
 pq.pop();

 a[source] -= n + 1;
 b[source] -= n + 1;

 if (a[source] > 0)
 pq.push(sample_pair(S(source)[a[source] - 1], source));
 }
 }
 }

 // postconditions: a[i] == b[i] in most cases, except when a[i] has been
 // clamped because of having reached the boundary

 // now return the result, calculate the offset, compare the keys on both
 // edges of the border

 // maximum of left edge, minimum of right edge
 ValueType* maxleft = NULL, * minright = NULL;
 for (diff_type i = 0; i < m; ++i)
 {
 if (a[i] > 0)
 {
 if (!maxleft)
 {
 maxleft = &(S(i)[a[i] - 1]);
 }
 else
 {
 // max
 if (comp(*maxleft, S(i)[a[i] - 1]))
 maxleft = &(S(i)[a[i] - 1]);
 }
 }
 if (b[i] < seqlen[i])
 {
 if (!minright)
 {
 minright = &(S(i)[b[i]]);
 }
 else
 {
 // min
 if (comp(S(i)[b[i]], *minright))
 minright = &(S(i)[b[i]]);
 }
 }
 }

 // minright is the splitter, in any case

 if (!maxleft || comp(*minright, *maxleft))
 {
 // good luck, everything is split unambigiously
 offset = 0;
 }
 else
 {
 // we have to calculate an offset
 offset = 0;

 for (diff_type i = 0; i < m; ++i)
 {
 diff_type lb = std::lower_bound(S(i), S(i) + seqlen[i], *minright, comp) - S(i);
 offset += a[i] - lb;
 }
 }

 delete[] seqlen;
 delete[] a;
 delete[] b;

 return *minright;
}

#undef S

} // namespace parallel

STXXL_END_NAMESPACE

#endif // !STXXL_PARALLEL_MULTISEQ_SELECTION_HEADER