637 lines
21 KiB
C++
637 lines
21 KiB
C++
|
/* boost random/discrete_distribution.hpp header file
|
||
|
|
*
|
||
|
|
* Copyright Steven Watanabe 2009-2011
|
||
|
|
* 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)
|
||
|
|
*
|
||
|
|
* See http://www.boost.org for most recent version including documentation.
|
||
|
|
*
|
||
|
|
* $Id$
|
||
|
|
*/
|
||
|
|
|
||
|
|
#ifndef BOOST_RANDOM_DISCRETE_DISTRIBUTION_HPP_INCLUDED
|
||
|
|
#define BOOST_RANDOM_DISCRETE_DISTRIBUTION_HPP_INCLUDED
|
||
|
|
|
||
|
|
#include <vector>
|
||
|
|
#include <limits>
|
||
|
|
#include <numeric>
|
||
|
|
#include <utility>
|
||
|
|
#include <iterator>
|
||
|
|
#include <boost/assert.hpp>
|
||
|
|
#include <boost/random/uniform_01.hpp>
|
||
|
|
#include <boost/random/uniform_int_distribution.hpp>
|
||
|
|
#include <boost/random/detail/config.hpp>
|
||
|
|
#include <boost/random/detail/operators.hpp>
|
||
|
|
#include <boost/random/detail/vector_io.hpp>
|
||
|
|
|
||
|
|
#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
|
||
|
|
#include <initializer_list>
|
||
|
|
#endif
|
||
|
|
|
||
|
|
#include <boost/range/begin.hpp>
|
||
|
|
#include <boost/range/end.hpp>
|
||
|
|
|
||
|
|
#include <boost/random/detail/disable_warnings.hpp>
|
||
|
|
|
||
|
|
namespace boost {
|
||
|
|
namespace random {
|
||
|
|
namespace detail {
|
||
|
|
|
||
|
|
template<class IntType, class WeightType>
|
||
|
|
struct integer_alias_table {
|
||
|
|
WeightType get_weight(IntType bin) const {
|
||
|
|
WeightType result = _average;
|
||
|
|
if(bin < _excess) ++result;
|
||
|
|
return result;
|
||
|
|
}
|
||
|
|
template<class Iter>
|
||
|
|
WeightType init_average(Iter begin, Iter end) {
|
||
|
|
WeightType weight_average = 0;
|
||
|
|
IntType excess = 0;
|
||
|
|
IntType n = 0;
|
||
|
|
// weight_average * n + excess == current partial sum
|
||
|
|
// This is a bit messy, but it's guaranteed not to overflow
|
||
|
|
for(Iter iter = begin; iter != end; ++iter) {
|
||
|
|
++n;
|
||
|
|
if(*iter < weight_average) {
|
||
|
|
WeightType diff = weight_average - *iter;
|
||
|
|
weight_average -= diff / n;
|
||
|
|
if(diff % n > excess) {
|
||
|
|
--weight_average;
|
||
|
|
excess += n - diff % n;
|
||
|
|
} else {
|
||
|
|
excess -= diff % n;
|
||
|
|
}
|
||
|
|
} else {
|
||
|
|
WeightType diff = *iter - weight_average;
|
||
|
|
weight_average += diff / n;
|
||
|
|
if(diff % n < n - excess) {
|
||
|
|
excess += diff % n;
|
||
|
|
} else {
|
||
|
|
++weight_average;
|
||
|
|
excess -= n - diff % n;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
_alias_table.resize(static_cast<std::size_t>(n));
|
||
|
|
_average = weight_average;
|
||
|
|
_excess = excess;
|
||
|
|
return weight_average;
|
||
|
|
}
|
||
|
|
void init_empty()
|
||
|
|
{
|
||
|
|
_alias_table.clear();
|
||
|
|
_alias_table.push_back(std::make_pair(static_cast<WeightType>(1),
|
||
|
|
static_cast<IntType>(0)));
|
||
|
|
_average = static_cast<WeightType>(1);
|
||
|
|
_excess = static_cast<IntType>(0);
|
||
|
|
}
|
||
|
|
bool operator==(const integer_alias_table& other) const
|
||
|
|
{
|
||
|
|
return _alias_table == other._alias_table &&
|
||
|
|
_average == other._average && _excess == other._excess;
|
||
|
|
}
|
||
|
|
static WeightType normalize(WeightType val, WeightType average)
|
||
|
|
{
|
||
|
|
return val;
|
||
|
|
}
|
||
|
|
static void normalize(std::vector<WeightType>&) {}
|
||
|
|
template<class URNG>
|
||
|
|
WeightType test(URNG &urng) const
|
||
|
|
{
|
||
|
|
return uniform_int_distribution<WeightType>(0, _average)(urng);
|
||
|
|
}
|
||
|
|
bool accept(IntType result, WeightType val) const
|
||
|
|
{
|
||
|
|
return result < _excess || val < _average;
|
||
|
|
}
|
||
|
|
static WeightType try_get_sum(const std::vector<WeightType>& weights)
|
||
|
|
{
|
||
|
|
WeightType result = static_cast<WeightType>(0);
|
||
|
|
for(typename std::vector<WeightType>::const_iterator
|
||
|
|
iter = weights.begin(), end = weights.end();
|
||
|
|
iter != end; ++iter)
|
||
|
|
{
|
||
|
|
if((std::numeric_limits<WeightType>::max)() - result > *iter) {
|
||
|
|
return static_cast<WeightType>(0);
|
||
|
|
}
|
||
|
|
result += *iter;
|
||
|
|
}
|
||
|
|
return result;
|
||
|
|
}
|
||
|
|
template<class URNG>
|
||
|
|
static WeightType generate_in_range(URNG &urng, WeightType max)
|
||
|
|
{
|
||
|
|
return uniform_int_distribution<WeightType>(
|
||
|
|
static_cast<WeightType>(0), max-1)(urng);
|
||
|
|
}
|
||
|
|
typedef std::vector<std::pair<WeightType, IntType> > alias_table_t;
|
||
|
|
alias_table_t _alias_table;
|
||
|
|
WeightType _average;
|
||
|
|
IntType _excess;
|
||
|
|
};
|
||
|
|
|
||
|
|
template<class IntType, class WeightType>
|
||
|
|
struct real_alias_table {
|
||
|
|
WeightType get_weight(IntType) const
|
||
|
|
{
|
||
|
|
return WeightType(1.0);
|
||
|
|
}
|
||
|
|
template<class Iter>
|
||
|
|
WeightType init_average(Iter first, Iter last)
|
||
|
|
{
|
||
|
|
std::size_t size = std::distance(first, last);
|
||
|
|
WeightType weight_sum =
|
||
|
|
std::accumulate(first, last, static_cast<WeightType>(0));
|
||
|
|
_alias_table.resize(size);
|
||
|
|
return weight_sum / size;
|
||
|
|
}
|
||
|
|
void init_empty()
|
||
|
|
{
|
||
|
|
_alias_table.clear();
|
||
|
|
_alias_table.push_back(std::make_pair(static_cast<WeightType>(1),
|
||
|
|
static_cast<IntType>(0)));
|
||
|
|
}
|
||
|
|
bool operator==(const real_alias_table& other) const
|
||
|
|
{
|
||
|
|
return _alias_table == other._alias_table;
|
||
|
|
}
|
||
|
|
static WeightType normalize(WeightType val, WeightType average)
|
||
|
|
{
|
||
|
|
return val / average;
|
||
|
|
}
|
||
|
|
static void normalize(std::vector<WeightType>& weights)
|
||
|
|
{
|
||
|
|
WeightType sum =
|
||
|
|
std::accumulate(weights.begin(), weights.end(),
|
||
|
|
static_cast<WeightType>(0));
|
||
|
|
for(typename std::vector<WeightType>::iterator
|
||
|
|
iter = weights.begin(),
|
||
|
|
end = weights.end();
|
||
|
|
iter != end; ++iter)
|
||
|
|
{
|
||
|
|
*iter /= sum;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
template<class URNG>
|
||
|
|
WeightType test(URNG &urng) const
|
||
|
|
{
|
||
|
|
return uniform_01<WeightType>()(urng);
|
||
|
|
}
|
||
|
|
bool accept(IntType, WeightType) const
|
||
|
|
{
|
||
|
|
return true;
|
||
|
|
}
|
||
|
|
static WeightType try_get_sum(const std::vector<WeightType>& weights)
|
||
|
|
{
|
||
|
|
return static_cast<WeightType>(1);
|
||
|
|
}
|
||
|
|
template<class URNG>
|
||
|
|
static WeightType generate_in_range(URNG &urng, WeightType)
|
||
|
|
{
|
||
|
|
return uniform_01<WeightType>()(urng);
|
||
|
|
}
|
||
|
|
typedef std::vector<std::pair<WeightType, IntType> > alias_table_t;
|
||
|
|
alias_table_t _alias_table;
|
||
|
|
};
|
||
|
|
|
||
|
|
template<bool IsIntegral>
|
||
|
|
struct select_alias_table;
|
||
|
|
|
||
|
|
template<>
|
||
|
|
struct select_alias_table<true> {
|
||
|
|
template<class IntType, class WeightType>
|
||
|
|
struct apply {
|
||
|
|
typedef integer_alias_table<IntType, WeightType> type;
|
||
|
|
};
|
||
|
|
};
|
||
|
|
|
||
|
|
template<>
|
||
|
|
struct select_alias_table<false> {
|
||
|
|
template<class IntType, class WeightType>
|
||
|
|
struct apply {
|
||
|
|
typedef real_alias_table<IntType, WeightType> type;
|
||
|
|
};
|
||
|
|
};
|
||
|
|
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* The class @c discrete_distribution models a \random_distribution.
|
||
|
|
* It produces integers in the range [0, n) with the probability
|
||
|
|
* of producing each value is specified by the parameters of the
|
||
|
|
* distribution.
|
||
|
|
*/
|
||
|
|
template<class IntType = int, class WeightType = double>
|
||
|
|
class discrete_distribution {
|
||
|
|
public:
|
||
|
|
typedef WeightType input_type;
|
||
|
|
typedef IntType result_type;
|
||
|
|
|
||
|
|
class param_type {
|
||
|
|
public:
|
||
|
|
|
||
|
|
typedef discrete_distribution distribution_type;
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Constructs a @c param_type object, representing a distribution
|
||
|
|
* with \f$p(0) = 1\f$ and \f$p(k|k>0) = 0\f$.
|
||
|
|
*/
|
||
|
|
param_type() : _probabilities(1, static_cast<WeightType>(1)) {}
|
||
|
|
/**
|
||
|
|
* If @c first == @c last, equivalent to the default constructor.
|
||
|
|
* Otherwise, the values of the range represent weights for the
|
||
|
|
* possible values of the distribution.
|
||
|
|
*/
|
||
|
|
template<class Iter>
|
||
|
|
param_type(Iter first, Iter last) : _probabilities(first, last)
|
||
|
|
{
|
||
|
|
normalize();
|
||
|
|
}
|
||
|
|
#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
|
||
|
|
/**
|
||
|
|
* If wl.size() == 0, equivalent to the default constructor.
|
||
|
|
* Otherwise, the values of the @c initializer_list represent
|
||
|
|
* weights for the possible values of the distribution.
|
||
|
|
*/
|
||
|
|
param_type(const std::initializer_list<WeightType>& wl)
|
||
|
|
: _probabilities(wl)
|
||
|
|
{
|
||
|
|
normalize();
|
||
|
|
}
|
||
|
|
#endif
|
||
|
|
/**
|
||
|
|
* If the range is empty, equivalent to the default constructor.
|
||
|
|
* Otherwise, the elements of the range represent
|
||
|
|
* weights for the possible values of the distribution.
|
||
|
|
*/
|
||
|
|
template<class Range>
|
||
|
|
explicit param_type(const Range& range)
|
||
|
|
: _probabilities(boost::begin(range), boost::end(range))
|
||
|
|
{
|
||
|
|
normalize();
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* If nw is zero, equivalent to the default constructor.
|
||
|
|
* Otherwise, the range of the distribution is [0, nw),
|
||
|
|
* and the weights are found by calling fw with values
|
||
|
|
* evenly distributed between \f$\mbox{xmin} + \delta/2\f$ and
|
||
|
|
* \f$\mbox{xmax} - \delta/2\f$, where
|
||
|
|
* \f$\delta = (\mbox{xmax} - \mbox{xmin})/\mbox{nw}\f$.
|
||
|
|
*/
|
||
|
|
template<class Func>
|
||
|
|
param_type(std::size_t nw, double xmin, double xmax, Func fw)
|
||
|
|
{
|
||
|
|
std::size_t n = (nw == 0) ? 1 : nw;
|
||
|
|
double delta = (xmax - xmin) / n;
|
||
|
|
BOOST_ASSERT(delta > 0);
|
||
|
|
for(std::size_t k = 0; k < n; ++k) {
|
||
|
|
_probabilities.push_back(fw(xmin + k*delta + delta/2));
|
||
|
|
}
|
||
|
|
normalize();
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Returns a vector containing the probabilities of each possible
|
||
|
|
* value of the distribution.
|
||
|
|
*/
|
||
|
|
std::vector<WeightType> probabilities() const
|
||
|
|
{
|
||
|
|
return _probabilities;
|
||
|
|
}
|
||
|
|
|
||
|
|
/** Writes the parameters to a @c std::ostream. */
|
||
|
|
BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
|
||
|
|
{
|
||
|
|
detail::print_vector(os, parm._probabilities);
|
||
|
|
return os;
|
||
|
|
}
|
||
|
|
|
||
|
|
/** Reads the parameters from a @c std::istream. */
|
||
|
|
BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
|
||
|
|
{
|
||
|
|
std::vector<WeightType> temp;
|
||
|
|
detail::read_vector(is, temp);
|
||
|
|
if(is) {
|
||
|
|
parm._probabilities.swap(temp);
|
||
|
|
}
|
||
|
|
return is;
|
||
|
|
}
|
||
|
|
|
||
|
|
/** Returns true if the two sets of parameters are the same. */
|
||
|
|
BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
|
||
|
|
{
|
||
|
|
return lhs._probabilities == rhs._probabilities;
|
||
|
|
}
|
||
|
|
/** Returns true if the two sets of parameters are different. */
|
||
|
|
BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
|
||
|
|
private:
|
||
|
|
/// @cond show_private
|
||
|
|
friend class discrete_distribution;
|
||
|
|
explicit param_type(const discrete_distribution& dist)
|
||
|
|
: _probabilities(dist.probabilities())
|
||
|
|
{}
|
||
|
|
void normalize()
|
||
|
|
{
|
||
|
|
impl_type::normalize(_probabilities);
|
||
|
|
}
|
||
|
|
std::vector<WeightType> _probabilities;
|
||
|
|
/// @endcond
|
||
|
|
};
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Creates a new @c discrete_distribution object that has
|
||
|
|
* \f$p(0) = 1\f$ and \f$p(i|i>0) = 0\f$.
|
||
|
|
*/
|
||
|
|
discrete_distribution()
|
||
|
|
{
|
||
|
|
_impl.init_empty();
|
||
|
|
}
|
||
|
|
/**
|
||
|
|
* Constructs a discrete_distribution from an iterator range.
|
||
|
|
* If @c first == @c last, equivalent to the default constructor.
|
||
|
|
* Otherwise, the values of the range represent weights for the
|
||
|
|
* possible values of the distribution.
|
||
|
|
*/
|
||
|
|
template<class Iter>
|
||
|
|
discrete_distribution(Iter first, Iter last)
|
||
|
|
{
|
||
|
|
init(first, last);
|
||
|
|
}
|
||
|
|
#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
|
||
|
|
/**
|
||
|
|
* Constructs a @c discrete_distribution from a @c std::initializer_list.
|
||
|
|
* If the @c initializer_list is empty, equivalent to the default
|
||
|
|
* constructor. Otherwise, the values of the @c initializer_list
|
||
|
|
* represent weights for the possible values of the distribution.
|
||
|
|
* For example, given the distribution
|
||
|
|
*
|
||
|
|
* @code
|
||
|
|
* discrete_distribution<> dist{1, 4, 5};
|
||
|
|
* @endcode
|
||
|
|
*
|
||
|
|
* The probability of a 0 is 1/10, the probability of a 1 is 2/5,
|
||
|
|
* the probability of a 2 is 1/2, and no other values are possible.
|
||
|
|
*/
|
||
|
|
discrete_distribution(std::initializer_list<WeightType> wl)
|
||
|
|
{
|
||
|
|
init(wl.begin(), wl.end());
|
||
|
|
}
|
||
|
|
#endif
|
||
|
|
/**
|
||
|
|
* Constructs a discrete_distribution from a Boost.Range range.
|
||
|
|
* If the range is empty, equivalent to the default constructor.
|
||
|
|
* Otherwise, the values of the range represent weights for the
|
||
|
|
* possible values of the distribution.
|
||
|
|
*/
|
||
|
|
template<class Range>
|
||
|
|
explicit discrete_distribution(const Range& range)
|
||
|
|
{
|
||
|
|
init(boost::begin(range), boost::end(range));
|
||
|
|
}
|
||
|
|
/**
|
||
|
|
* Constructs a discrete_distribution that approximates a function.
|
||
|
|
* If nw is zero, equivalent to the default constructor.
|
||
|
|
* Otherwise, the range of the distribution is [0, nw),
|
||
|
|
* and the weights are found by calling fw with values
|
||
|
|
* evenly distributed between \f$\mbox{xmin} + \delta/2\f$ and
|
||
|
|
* \f$\mbox{xmax} - \delta/2\f$, where
|
||
|
|
* \f$\delta = (\mbox{xmax} - \mbox{xmin})/\mbox{nw}\f$.
|
||
|
|
*/
|
||
|
|
template<class Func>
|
||
|
|
discrete_distribution(std::size_t nw, double xmin, double xmax, Func fw)
|
||
|
|
{
|
||
|
|
std::size_t n = (nw == 0) ? 1 : nw;
|
||
|
|
double delta = (xmax - xmin) / n;
|
||
|
|
BOOST_ASSERT(delta > 0);
|
||
|
|
std::vector<WeightType> weights;
|
||
|
|
for(std::size_t k = 0; k < n; ++k) {
|
||
|
|
weights.push_back(fw(xmin + k*delta + delta/2));
|
||
|
|
}
|
||
|
|
init(weights.begin(), weights.end());
|
||
|
|
}
|
||
|
|
/**
|
||
|
|
* Constructs a discrete_distribution from its parameters.
|
||
|
|
*/
|
||
|
|
explicit discrete_distribution(const param_type& parm)
|
||
|
|
{
|
||
|
|
param(parm);
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Returns a value distributed according to the parameters of the
|
||
|
|
* discrete_distribution.
|
||
|
|
*/
|
||
|
|
template<class URNG>
|
||
|
|
IntType operator()(URNG& urng) const
|
||
|
|
{
|
||
|
|
BOOST_ASSERT(!_impl._alias_table.empty());
|
||
|
|
IntType result;
|
||
|
|
WeightType test;
|
||
|
|
do {
|
||
|
|
result = uniform_int_distribution<IntType>((min)(), (max)())(urng);
|
||
|
|
test = _impl.test(urng);
|
||
|
|
} while(!_impl.accept(result, test));
|
||
|
|
if(test < _impl._alias_table[static_cast<std::size_t>(result)].first) {
|
||
|
|
return result;
|
||
|
|
} else {
|
||
|
|
return(_impl._alias_table[static_cast<std::size_t>(result)].second);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Returns a value distributed according to the parameters
|
||
|
|
* specified by param.
|
||
|
|
*/
|
||
|
|
template<class URNG>
|
||
|
|
IntType operator()(URNG& urng, const param_type& parm) const
|
||
|
|
{
|
||
|
|
if(WeightType limit = impl_type::try_get_sum(parm._probabilities)) {
|
||
|
|
WeightType val = impl_type::generate_in_range(urng, limit);
|
||
|
|
WeightType sum = 0;
|
||
|
|
std::size_t result = 0;
|
||
|
|
for(typename std::vector<WeightType>::const_iterator
|
||
|
|
iter = parm._probabilities.begin(),
|
||
|
|
end = parm._probabilities.end();
|
||
|
|
iter != end; ++iter, ++result)
|
||
|
|
{
|
||
|
|
sum += *iter;
|
||
|
|
if(sum > val) {
|
||
|
|
return result;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
// This shouldn't be reachable, but round-off error
|
||
|
|
// can prevent any match from being found when val is
|
||
|
|
// very close to 1.
|
||
|
|
return static_cast<IntType>(parm._probabilities.size() - 1);
|
||
|
|
} else {
|
||
|
|
// WeightType is integral and sum(parm._probabilities)
|
||
|
|
// would overflow. Just use the easy solution.
|
||
|
|
return discrete_distribution(parm)(urng);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
/** Returns the smallest value that the distribution can produce. */
|
||
|
|
result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
|
||
|
|
/** Returns the largest value that the distribution can produce. */
|
||
|
|
result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
|
||
|
|
{ return static_cast<result_type>(_impl._alias_table.size() - 1); }
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Returns a vector containing the probabilities of each
|
||
|
|
* value of the distribution. For example, given
|
||
|
|
*
|
||
|
|
* @code
|
||
|
|
* discrete_distribution<> dist = { 1, 4, 5 };
|
||
|
|
* std::vector<double> p = dist.param();
|
||
|
|
* @endcode
|
||
|
|
*
|
||
|
|
* the vector, p will contain {0.1, 0.4, 0.5}.
|
||
|
|
*
|
||
|
|
* If @c WeightType is integral, then the weights
|
||
|
|
* will be returned unchanged.
|
||
|
|
*/
|
||
|
|
std::vector<WeightType> probabilities() const
|
||
|
|
{
|
||
|
|
std::vector<WeightType> result(_impl._alias_table.size(), static_cast<WeightType>(0));
|
||
|
|
std::size_t i = 0;
|
||
|
|
for(typename impl_type::alias_table_t::const_iterator
|
||
|
|
iter = _impl._alias_table.begin(),
|
||
|
|
end = _impl._alias_table.end();
|
||
|
|
iter != end; ++iter, ++i)
|
||
|
|
{
|
||
|
|
WeightType val = iter->first;
|
||
|
|
result[i] += val;
|
||
|
|
result[static_cast<std::size_t>(iter->second)] += _impl.get_weight(i) - val;
|
||
|
|
}
|
||
|
|
impl_type::normalize(result);
|
||
|
|
return(result);
|
||
|
|
}
|
||
|
|
|
||
|
|
/** Returns the parameters of the distribution. */
|
||
|
|
param_type param() const
|
||
|
|
{
|
||
|
|
return param_type(*this);
|
||
|
|
}
|
||
|
|
/** Sets the parameters of the distribution. */
|
||
|
|
void param(const param_type& parm)
|
||
|
|
{
|
||
|
|
init(parm._probabilities.begin(), parm._probabilities.end());
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Effects: Subsequent uses of the distribution do not depend
|
||
|
|
* on values produced by any engine prior to invoking reset.
|
||
|
|
*/
|
||
|
|
void reset() {}
|
||
|
|
|
||
|
|
/** Writes a distribution to a @c std::ostream. */
|
||
|
|
BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, discrete_distribution, dd)
|
||
|
|
{
|
||
|
|
os << dd.param();
|
||
|
|
return os;
|
||
|
|
}
|
||
|
|
|
||
|
|
/** Reads a distribution from a @c std::istream */
|
||
|
|
BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, discrete_distribution, dd)
|
||
|
|
{
|
||
|
|
param_type parm;
|
||
|
|
if(is >> parm) {
|
||
|
|
dd.param(parm);
|
||
|
|
}
|
||
|
|
return is;
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Returns true if the two distributions will return the
|
||
|
|
* same sequence of values, when passed equal generators.
|
||
|
|
*/
|
||
|
|
BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(discrete_distribution, lhs, rhs)
|
||
|
|
{
|
||
|
|
return lhs._impl == rhs._impl;
|
||
|
|
}
|
||
|
|
/**
|
||
|
|
* Returns true if the two distributions may return different
|
||
|
|
* sequences of values, when passed equal generators.
|
||
|
|
*/
|
||
|
|
BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(discrete_distribution)
|
||
|
|
|
||
|
|
private:
|
||
|
|
|
||
|
|
/// @cond show_private
|
||
|
|
|
||
|
|
template<class Iter>
|
||
|
|
void init(Iter first, Iter last, std::input_iterator_tag)
|
||
|
|
{
|
||
|
|
std::vector<WeightType> temp(first, last);
|
||
|
|
init(temp.begin(), temp.end());
|
||
|
|
}
|
||
|
|
template<class Iter>
|
||
|
|
void init(Iter first, Iter last, std::forward_iterator_tag)
|
||
|
|
{
|
||
|
|
std::vector<std::pair<WeightType, IntType> > below_average;
|
||
|
|
std::vector<std::pair<WeightType, IntType> > above_average;
|
||
|
|
WeightType weight_average = _impl.init_average(first, last);
|
||
|
|
WeightType normalized_average = _impl.get_weight(0);
|
||
|
|
std::size_t i = 0;
|
||
|
|
for(; first != last; ++first, ++i) {
|
||
|
|
WeightType val = impl_type::normalize(*first, weight_average);
|
||
|
|
std::pair<WeightType, IntType> elem(val, static_cast<IntType>(i));
|
||
|
|
if(val < normalized_average) {
|
||
|
|
below_average.push_back(elem);
|
||
|
|
} else {
|
||
|
|
above_average.push_back(elem);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
typename impl_type::alias_table_t::iterator
|
||
|
|
b_iter = below_average.begin(),
|
||
|
|
b_end = below_average.end(),
|
||
|
|
a_iter = above_average.begin(),
|
||
|
|
a_end = above_average.end()
|
||
|
|
;
|
||
|
|
while(b_iter != b_end && a_iter != a_end) {
|
||
|
|
_impl._alias_table[static_cast<std::size_t>(b_iter->second)] =
|
||
|
|
std::make_pair(b_iter->first, a_iter->second);
|
||
|
|
a_iter->first -= (_impl.get_weight(b_iter->second) - b_iter->first);
|
||
|
|
if(a_iter->first < normalized_average) {
|
||
|
|
*b_iter = *a_iter++;
|
||
|
|
} else {
|
||
|
|
++b_iter;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
for(; b_iter != b_end; ++b_iter) {
|
||
|
|
_impl._alias_table[static_cast<std::size_t>(b_iter->second)].first =
|
||
|
|
_impl.get_weight(b_iter->second);
|
||
|
|
}
|
||
|
|
for(; a_iter != a_end; ++a_iter) {
|
||
|
|
_impl._alias_table[static_cast<std::size_t>(a_iter->second)].first =
|
||
|
|
_impl.get_weight(a_iter->second);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
template<class Iter>
|
||
|
|
void init(Iter first, Iter last)
|
||
|
|
{
|
||
|
|
if(first == last) {
|
||
|
|
_impl.init_empty();
|
||
|
|
} else {
|
||
|
|
typename std::iterator_traits<Iter>::iterator_category category;
|
||
|
|
init(first, last, category);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
typedef typename detail::select_alias_table<
|
||
|
|
(::boost::is_integral<WeightType>::value)
|
||
|
|
>::template apply<IntType, WeightType>::type impl_type;
|
||
|
|
impl_type _impl;
|
||
|
|
/// @endcond
|
||
|
|
};
|
||
|
|
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
#include <boost/random/detail/enable_warnings.hpp>
|
||
|
|
|
||
|
|
#endif
|