Random Generation From Sequences

Question

With the inclusion of <random> into C++11, we can finally chuck out std::rand and start using generator with much better properties. Still, it's easy to get things wrong (uniform sampling where 0 occurs more than it should due to usage of % size comes to mind). To that end, I've written some wrapper libraries to ease usage of the library for some basic tasks, namely:

• Selection of values from a sequence with equal probability
• Selection of values from a sequence according to a weighted distribution
• Random "choice" - k out of n selection.

The style is very similar to that of most of the algorithms in <algorithm>. Any comments or suggestions (or bug-finding) is welcome:

initialized_generator.hpp:

/*! \file intialized_generator.hpp
*   \brief Helper struct to initialize and "warm up" a random number generator.
*
*/

#ifndef INITIALIZED_GENERATOR_SGL_HPP__
#define INITIALIZED_GENERATOR_SGL_HPP__

#include <random>
#include <array>
#include <algorithm>
#include <type_traits>
#include <ctime>

namespace simplegl
{
namespace detail
{
template <typename GeneratorType>
struct initialized_generator
{
public:

    typedef GeneratorType                        generator_type;
    typedef typename generator_type::result_type result_type;

private:

    generator_type                               generator_;

    //! \brief Seeds the given generator with a given source of randomness.
    /*!  
     *  Some random number generators generate "poor" random numbers until their
     *  internal states are sufficiently "mixed up". This should seed the given
     *  generator with enough randomness to overcome the initial bad statistical
     *  properties that are seen in some of these generators.
     *
     *  The code below was nabbed from 
     *  http://stackoverflow.com/questions/15509270/does-stdmt19937-require-warmup 
     *
     *  \param rd A std::random_device, utilized to generate random seed data.
    */
    void warmup(std::random_device& rd)
    {
        std::array<int, std::mt19937::state_size> seed_data;
        std::generate_n(seed_data.data(), seed_data.size(), std::ref(rd));
        std::seed_seq seq(std::begin(seed_data), std::end(seed_data));
        generator_.seed(seq);
    }

    //! \brief Seeds the given generator utilizing the current time.
    /*!  
     *  Seeds the given generator utilizing the current time as a seed.
     *  If the generator is std::mt19937, then also attempts to initialize
     *  it and warm it up with some additional (not really random) seed
     *  parameters.    
     */
    void warmup()
    {
        std::time_t initial = time(nullptr);
        //This branch will be optimized out at compile time
        if(std::is_same<GeneratorType, std::mt19937>::value) {
            //Even relatively close "non-random" numbers for a sequence should
            //be ok for seeding a mersenne twister.
            std::seed_seq seq{initial, initial + 1, initial + 2, initial + 3,
                              initial + 4, initial + 5};
            generator_.seed(seq);
        } else {
            //For an LCG, time should suffice as a seed. For a subtract with carry,
            //it may not - apparently they are actually quite difficult to seed
            //and are incredibly sensitive to initial state. This should perhaps
            //be changed to deal with that, but I don't really have the expertise
            //to know exactly how.
            generator_.seed(initial);
        }
    }

public:

    initialized_generator()
    {
        warmup();
    }

    initialized_generator(std::random_device& rd)
    {
        warmup(rd);
    }

    //! \brief Returns the minimum possible value the underlying generator 
    /*!        can generate.
     *
     *  Simple wrapper function forwarding to GeneratorType::min().   
     */
    static constexpr result_type min()
    {
        return generator_type::min();
    }

    //! \brief Returns the maximum possible value the underlying generator 
    /*!        can generate.
     *
     *  Simple wrapper function forwarding to GeneratorType::max().   
     */
    static constexpr result_type max()
    {
        return generator_type::max();
    }

    //! \brief Returns the next random number from the underlying generator.
    /*! 
     *  Simple wrapper function forwarding to operator()().
     *  \return The next random number, of type GeneratorType::result_type.
     */
    result_type operator()()
    {
        return generator_();
    }

    //! \brief Allows implicit conversion back to the underlying GeneratorType.
    /*!  
     */
    operator generator_type()
    {
        return generator_;
    }

}; //end struct initialized_generator

} //end namespace detail 
} //end namespace simplegl

#endif //INITIALIZED_GENERATOR_SGL_HPP__

---

weighted_selection.hpp:

/*! \file weighted_selection.hpp
*   \brief Selects values from a sequence based on a linear weighting.
*
*/

#ifndef WEIGHTED_SELECTION_SGL_HPP__
#define WEIGHTED_SELECTION_SGL_HPP__

#include <iterator>
#include <random>
#include <initializer_list>
#include <utility>

#include "initialized_generator.hpp"

namespace simplegl
{

template <typename GeneratorType = std::mt19937, typename UIntType = std::size_t>
struct weighted_selection
{
private:

    typedef UIntType uint_type;

    std::discrete_distribution<UIntType>         distribution_;
    detail::initialized_generator<GeneratorType> generator_;
    uint_type                                    prob_length_;

public:

    weighted_selection(std::initializer_list<double> init, std::random_device& rd)
      : distribution_(init),
        generator_(rd),
        prob_length_(init.size())
    { }

    weighted_selection(std::initializer_list<double> init)
      : distribution_(init),
        prob_length_(init.size())
    { }

    template <typename Iterator>
    weighted_selection(Iterator begin, Iterator end, std::random_device& rd)
      : distribution_(begin, end),
        generator_(rd),
        prob_length_(std::distance(begin, end))
    { }

    template <typename Iterator>
    weighted_selection(Iterator begin, Iterator end)
      : distribution_(begin, end),
        prob_length_(std::distance(begin, end))
    { }

    //! \brief Selects a random value from the sequence [begin, end) based
    //*        on the weights given on initialization.
    /*!  
     *  Overload that performs error checking based on the distance
     *  between begin and end.
     *
     *  \param begin A parameter supporting the RandomAccessIterator concept,
     *         which dereferences to the start of the given sequenece.
     *  \param end A parameter supporting the RandomAccessIterator concept,
     *         which is one past the end of the sequence.
     *  \return The next random value, with probability given by the 
     *          supplied weights. 
     *  \throw std::length_error when the distance between end and
     *         begin is not equal to the length of the weights initially supplied
     *         on object construction.
    */
    template <typename RandomAccessIter>
    typename std::iterator_traits<RandomAccessIter>::value_type
    operator()(RandomAccessIter begin, RandomAccessIter end)
    {
        auto distance = end - begin;
        uint_type unsigned_distance = 
            (distance > 0) ? static_cast<uint_type>(distance) : 0;
        if(unsigned_distance != prob_length_) {
            throw std::length_error("Iterator distance does not equal probability \
                                     weight distances");
        }
        auto next = distribution_(generator_);
        return *(begin + next);
    }

    //! \brief Selects a random value from the sequence [begin, ...) based
    //*        on the weights given on initialization.
    /*!  
     *  Overloaded version that performs no error checking. It is assumed that
     *  (begin + weights.size() - 1) is dereferencable.
     *
     *  \param begin A parameter supporting the RandomAccessIterator concept,
     *         which dereferences to the start of the given sequenece.
     *  \return The next random value, with probability given by the 
     *          supplied weights. 
    */
    template <typename RandomAccessIter>
    typename std::iterator_traits<RandomAccessIter>::value_type
    operator()(RandomAccessIter begin)
    {
        auto next = distribution_(generator_);
        return *(begin + next);
    }

    //! \brief The number of weights this was initialized with.
    uint_type num_weights() const
    {
        return prob_length_;
    }

}; //end struct weighted_selection

} //end namespace simplegl

#endif //WEIGHTED_SELECTION_SGL_HPP__

---

uniform_selection.hpp:

/*! \file uniform_selection.hpp
*   \brief Selects values from a sequence with equal probability.
*
*/

#ifndef UNIFORM_SELECTION_SGL_HPP__
#define UNIFORM_SELECTION_SGL_HPP__

#include <iterator>
#include <random>
#include <stdexcept>

#include "initialized_generator.hpp"

namespace simplegl
{

template <typename GeneratorType = std::mt19937, typename UIntType = std::size_t>
struct uniform_selection
{
public:

    typedef UIntType uint_type;

private:

    std::uniform_int_distribution<UIntType>      distribution_;
    detail::initialized_generator<GeneratorType> generator_;

public:

    uniform_selection()
    { }

    uniform_selection(std::random_device& rd)
      : generator_(rd)
    { }

    //! \brief Selects a random value from the sequence [begin, end) with
    //*        equal probability.
    /*!  
     *  \param begin A parameter supporting the RandomAccessIterator concept,
     *         which dereferences to the start of the given sequenece.
     *  \param end A parameter suporting the RandomAccessIterator concept,
     *         which points to one past the end of the sequence.
     *  \return The next random value from [begin, end) with each value
     *          having probability of 1 / distance(begin, end) of being
     *          chosen.
     *  \throw  std::length_error if begin and end are equal, or if
     *          the distance between them is negative.
    */     
    template <typename RandomAccessIter>
    typename std::iterator_traits<RandomAccessIter>::value_type
    operator()(RandomAccessIter begin, RandomAccessIter end)
    {
        typedef typename std::uniform_int_distribution<UIntType> distribution_type;
        typedef typename distribution_type::param_type param_type;

        auto distance = end - begin;
        uint_type unsigned_distance = 
            (distance > 0) ? static_cast<uint_type>(distance) : 0;
        uint_type zero = 0;

        if(distance <= 0) {
            throw std::length_error("Iterator distance is negative or zero.");
        }

        auto next = distribution_(generator_, param_type(zero, unsigned_distance - 1));
        return *(begin + next);
    }
}; //end struct uniform_selection

} //end namespace simplegl

#endif //UNIFORM_SELECTION_SGL_HPP__

---

intrusive_random_choice.hpp:

/*! \file intrusive_random_choice.hpp
*   \brief Selects k random values from a sequence of length n, k <= n.
*
*/

#ifndef INTRUSIVE_RANDOM_CHOICE_SGL_HPP__
#define INTRUSIVE_RANDOM_CHOICE_SGL_HPP__

#include <algorithm>
#include <random>
#include <limits>

#include "initialized_generator.hpp"

namespace simplegl
{

template <typename GeneratorType = std::mt19937, typename UIntType = std::size_t>
struct intrusive_random_choice
{
public:

    typedef UIntType uint_type;

private:

    detail::initialized_generator<GeneratorType> generator_;

public:

    intrusive_random_choice()
    { }

    intrusive_random_choice(std::random_device& rd)
      : generator_(rd)
    { }

    //! \brief Selects k random values from a sequence, where each value is 
    //*        chosen with equal probability. This will modify the order in the
    //*        underlying sequence [begin, end).
    /*!
     *  k must be less than or equal to the length of the sequence. No checking
     *  is done to enforce this. If k is greater than the sequence length, this
     *  will result in undefined behaviour. Inserts the k values into a sequence
     *  starting at insert. There must be space allocated for at least k elements,
     *  otherwise, undefined behaviour will result. 
     *
     *  This will generate a correct uniform distribution for k when the sequence 
     *  length is less than (roughly) 2500. Any size larger requires more
     *  sophisticated techniques.
     *     
     *  \param begin A parameter supporting the RandomAccessIterator concept,
     *         which dereferences to the start of the given sequenece to select
     *         values from.
     *  \param end A parameter suporting the RandomAccessIterator concept,
     *         which points to one past the end of the sequence to select
     *         values from.
     *  \param insert The start of the sequence where the k random values will
     *         be inserted. 
    */  
    template <typename RandomAccessIter, typename RandomAccessIter2>
    void operator()(RandomAccessIter begin, RandomAccessIter end, 
                    RandomAccessIter2 insert, uint_type k)
    {
        std::shuffle(begin, end, generator_);
        for(uint_type i = 0; i < k; ++i) {
            *insert++ = *begin++;
        }
    }

}; //end struct intrusive_random_choice

} //end namespace simplegl

#endif //INTRUSIVE_RANDOM_CHOICE_SGL_HPP__

---

example_usage.cpp:

#include <initializer_list>
#include <vector>
#include <string>
#include <iostream>
#include <unordered_map>
#include <deque>

#include "weighted_selection.hpp"
#include "uniform_selection.hpp"
#include "intrusive_random_choice.hpp"

using namespace simplegl;

int main()
{
    //Weighted Selection example
    weighted_selection<> rs {10, 15, 20, 12, 18};

    std::vector<std::string> s {"The", "Cake", "Is", "A", "Lie"};
    std::unordered_map<std::string, std::size_t> strcount;

    //Select 15 random elements:
    for(int i = 0; i < 15; ++i) {
        std::cout << rs(s.begin(), s.end()) << "\n";
    }

    std::cout << "\n--------------------------------------------------\n";

    //Uniform Selection Example
    uniform_selection<> us;

    std::deque<int> d{10, 15, 20, 25, 30};
    std::unordered_map<int, std::size_t> counts;

    for(unsigned i = 0; i < 100000; ++i) {
        ++counts[us(d.begin(), d.end())];
    }

    for(const auto& it : counts) {
        std::cout << it.first << " seen: " << it.second << " times\n";
    }

    std::cout << "\n";

    //The same uniform_selection object can be used with sequences with
    //different lengths
    std::vector<int> ex {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
    for(int i = 0; i < 5; ++i) {
        std::cout << "Random element from ex: " << us(ex.begin(), ex.end()) << "\n";
    }

    std::cout << "\n--------------------------------------------------\n";

    //Intrusive random choice example
    intrusive_random_choice<> irc;
    //Make a copy of d so we leave it in order
    std::deque<int> dcopy(d.begin(), d.end());
    //We'll select 3 elements from the 5 in dcopy
    std::vector<int> v(3);
    counts.clear();

    for(int j = 0; j < 100000; ++ j) {
        irc(dcopy.begin(), dcopy.end(), v.begin(), 3);
        for(int i : v) {
            ++counts[i];
        }
    }

    for(const auto& it : counts) {
        std::cout << it.first << " seen: " << it.second << " times\n";
    }
}

Answer 1

Overall, your code is clean, well written and seems very consistent (and the comments also seem relevant). There is not much to say, so I will focus on some details since I don't see anything else to review:

• First of all, Loki is right in his comment: your header guards names (UNIFORM_SELECTION_SGL_HPP__, INTRUSIVE_RANDOM_CHOICE_SGL_HPP__...) are actually reserved to the implementation since they contain double underscores. You should remove the last underscore so that they are well-formed.
• It is mainly a matter of taste, but I would use the new type alias syntax (with using) instead of the old typedef. Contrary to typedef, using can be templated, therefore it allows you to always be consistent, whether you have templates or not. Also, I find the X = Y syntax clearer and closer to variable assignment.

• A tidbit:

  std::is_same<GeneratorType, std::mt19937>::value
  

  I would replace GeneratorType by generator_type. It won't change the generated code but since you use the typedefs instead of the template parameters names in the rest of your code, it will be more consistent.

• //This branch will be optimized out at compile time

  I would replace will by should. While any decent compiler should optimize it, there is no such garantee that it will be the case for any compiler. Since you don't know which compiler your code will be compiled with, you shouldn't make the assumption (Ok, that's really a detail).

• You could use curly braces instead of parenthesis in the initialization list of your constructors to prevent implicit type conversions.

As you can see, there is almost nothing to say about your current code. You even used std::shuffle and that's great (even greater since std::random_shuffle will be deprecated in C++14). That's some really good code, kudos! :)