[Boost-commit] svn:boost r63180 - in trunk: boost/random libs/random/doc libs/random/test

Subject: [Boost-commit] svn:boost r63180 - in trunk: boost/random libs/random/doc libs/random/test
From: steven_at_[hidden]
Date: 2010-06-21 00:03:57

Author: steven_watanabe
Date: 2010-06-21 00:03:55 EDT (Mon, 21 Jun 2010)
New Revision: 63180
URL: http://svn.boost.org/trac/boost/changeset/63180

Log:
Implement discrete_distribution.  Fixes #920.
Added:
   trunk/boost/random/discrete_distribution.hpp   (contents, props changed)
   trunk/libs/random/test/test_discrete.cpp   (contents, props changed)
   trunk/libs/random/test/test_discrete_distribution.cpp   (contents, props changed)
Text files modified: 
   trunk/libs/random/doc/Jamfile.v2        |     1 +                                       
   trunk/libs/random/doc/distributions.qbk |     1 +                                       
   trunk/libs/random/doc/random.qbk        |     1 +                                       
   trunk/libs/random/test/Jamfile.v2       |     2 ++                                      
   4 files changed, 5 insertions(+), 0 deletions(-)

Added: trunk/boost/random/discrete_distribution.hpp
==============================================================================
--- (empty file)
+++ trunk/boost/random/discrete_distribution.hpp	2010-06-21 00:03:55 EDT (Mon, 21 Jun 2010)
@@ -0,0 +1,494 @@
+// discrete_distribution.hpp
+//
+// Copyright (c) 2009-2010
+// Steven Watanabe
+//
+// 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 BOOST_RANDOM_DISCRETE_DISTRIBUTION_HPP_INCLUDED
+#define BOOST_RANDOM_DISCRETE_DISTRIBUTION_HPP_INCLUDED
+
+#include <vector>
+#include <cassert>
+#include <limits>
+#include <numeric>
+#include <utility>
+#include <iterator>
+#include <boost/random/uniform_01.hpp>
+#include <boost/random/uniform_int.hpp>
+#include <boost/random/detail/config.hpp>
+
+#ifndef BOOST_NO_INITIALIZER_LISTS
+#include <initializer_list>
+#endif
+
+#include <boost/range/begin.hpp>
+#include <boost/range/end.hpp>
+
+namespace boost {
+namespace random {
+
+/**
+ * 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:
+        /**
+         * 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_INITIALIZER_LISTS
+        /**
+         * 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;
+            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;
+        }
+
+#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
+        /** Writes the parameters to a @c std::ostream. */
+        template<class CharT, class Traits>
+        friend std::basic_ostream<CharT, Traits>&
+        operator<<(std::basic_ostream<CharT, Traits>& os,
+                   const param_type& parm)
+        {
+            typename std::vector<WeightType>::const_iterator
+                iter = parm._probabilities.begin(),
+                end =  parm._probabilities.end();
+            os << '[';
+            if(iter != end) {
+                os << *iter;
+                ++iter;
+                for(; iter != end; ++iter)
+                {
+                    os << ' ' << *iter;
+                }
+            }
+            os << ']';
+            return os;
+        }
+        
+        /** Reads the parameters from a @c std::istream. */
+        template<class CharT, class Traits>
+        friend std::basic_istream<CharT, Traits>&
+        operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm)
+        {
+            std::vector<WeightType> result;
+            char ch;
+            if(!(is >> ch)) {
+                return is;
+            }
+            if(ch != '[') {
+                is.putback(ch);
+                is.setstate(std::ios_base::failbit);
+                return is;
+            }
+            WeightType val;
+            while(is >> std::ws >> val) {
+                result.push_back(val);
+            }
+            if(is.fail()) {
+                is.clear();
+                if(!(is >> ch)) {
+                    return is;
+                }
+                if(ch != ']') {
+                    is.putback(ch);
+                    is.setstate(std::ios_base::failbit);
+                }
+            }
+            parm._probabilities.swap(result);
+            return is;
+        }
+#endif
+
+        /** Returns true if the two sets of parameters are the same. */
+        friend bool operator==(const param_type& lhs, const param_type& rhs)
+        {
+            return lhs._probabilities == rhs._probabilities;
+        }
+        /** Returns true if the two sets of parameters are different. */
+        friend bool operator!=(const param_type& lhs, const param_type& rhs)
+        {
+            return !(lhs == rhs);
+        }
+    private:
+        /// @cond
+        friend class discrete_distribution;
+        explicit param_type(const discrete_distribution& dist)
+          : _probabilities(dist.probabilities())
+        {}
+        void normalize()
+        {
+            WeightType sum =
+                std::accumulate(_probabilities.begin(), _probabilities.end(),
+                                static_cast<WeightType>(0));
+            for(typename std::vector<WeightType>::iterator
+                    iter = _probabilities.begin(),
+                    end = _probabilities.end();
+                    iter != end; ++iter)
+            {
+                *iter /= sum;
+            }
+        }
+        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()
+    {
+        _alias_table.push_back(std::make_pair(static_cast<WeightType>(1),
+                                              static_cast<IntType>(0)));
+    }
+    /**
+     * 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_INITIALIZER_LISTS
+    /**
+     * 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;
+        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 {
+        assert(!_alias_table.empty());
+        WeightType test = uniform_01<WeightType>()(urng);
+        IntType result = uniform_int<IntType>((min)(), (max)())(urng);
+        if(test < _alias_table[result].first) {
+            return result;
+        } else {
+            return(_alias_table[result].second);
+        }
+    }
+    
+    /**
+     * Returns a value distributed according to the parameters
+     * specified by parm.
+     */
+    template<class URNG>
+    IntType operator()(URNG& urng, const param_type& parm) const
+    {
+        while(true) {
+            WeightType val = uniform_01<WeightType>()(urng);
+            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;
+                }
+            }
+        }
+    }
+    
+    /** 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>(_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}.
+     */
+    std::vector<WeightType> probabilities() const
+    {
+        std::vector<WeightType> result(_alias_table.size());
+        const WeightType mean =
+            static_cast<WeightType>(1) / _alias_table.size();
+        std::size_t i = 0;
+        for(typename alias_table_t::const_iterator
+                iter = _alias_table.begin(),
+                end = _alias_table.end();
+                iter != end; ++iter, ++i)
+        {
+            WeightType val = iter->first * mean;
+            result[i] += val;
+            result[iter->second] += mean - val;
+        }
+        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() {}
+
+#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
+    /** Writes a distribution to a @c std::ostream. */
+    template<class CharT, class Traits>
+    friend std::basic_ostream<CharT, Traits>&
+    operator<<(std::basic_ostream<CharT, Traits>& os,
+               const discrete_distribution& dd)
+    {
+        os << dd.param();
+        return os;
+    }
+
+    /** Reads a distribution from a @c std::istream */
+    template<class CharT, class Traits>
+    friend std::basic_istream<CharT, Traits>&
+    operator>>(std::basic_istream<CharT, Traits>& is, discrete_distribution& dd)
+    {
+        param_type parm;
+        if(is >> parm) {
+            dd.param(parm);
+        }
+        return is;
+    }
+#endif
+
+    /**
+     * Returns true if the two distributions will return the
+     * same sequence of values, when passed equal generators.
+     */
+    friend bool operator==(const discrete_distribution& lhs,
+                           const discrete_distribution& rhs)
+    {
+        return lhs._alias_table == rhs._alias_table;
+    }
+    /**
+     * Returns true if the two distributions may return different
+     * sequences of values, when passed equal generators.
+     */
+    friend bool operator!=(const discrete_distribution& lhs,
+                           const discrete_distribution& rhs)
+    {
+        return !(lhs == rhs);
+    }
+
+private:
+
+    /// @cond
+
+    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;
+        std::size_t size = std::distance(first, last);
+        WeightType weight_sum =
+            std::accumulate(first, last, static_cast<WeightType>(0));
+        WeightType weight_average = weight_sum / size;
+        std::size_t i = 0;
+        for(; first != last; ++first, ++i) {
+            WeightType val = *first / weight_average;
+            std::pair<WeightType, IntType> elem(val, static_cast<IntType>(i));
+            if(val < static_cast<WeightType>(1)) {
+                below_average.push_back(elem);
+            } else {
+                above_average.push_back(elem);
+            }
+        }
+
+        _alias_table.resize(size);
+        typename 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) {
+            _alias_table[b_iter->second] =
+                std::make_pair(b_iter->first, a_iter->second);
+            a_iter->first -= (static_cast<WeightType>(1) - b_iter->first);
+            if(a_iter->first < static_cast<WeightType>(1)) {
+                *b_iter = *a_iter++;
+            } else {
+                ++b_iter;
+            }
+        }
+        for(; b_iter != b_end; ++b_iter) {
+            _alias_table[b_iter->second].first = static_cast<WeightType>(1);
+        }
+        for(; a_iter != a_end; ++a_iter) {
+            _alias_table[a_iter->second].first = static_cast<WeightType>(1);
+        }
+    }
+    template<class Iter>
+    void init(Iter first, Iter last)
+    {
+        if(first == last) {
+            _alias_table.clear();
+            _alias_table.push_back(std::make_pair(static_cast<WeightType>(1),
+                                                  static_cast<IntType>(0)));
+        } else {
+            typename std::iterator_traits<Iter>::iterator_category category;
+            init(first, last, category);
+        }
+    }
+    typedef std::vector<std::pair<WeightType, IntType> > alias_table_t;
+    alias_table_t _alias_table;
+    /// @endcond
+};
+
+}
+}
+
+#endif

Modified: trunk/libs/random/doc/Jamfile.v2
==============================================================================
--- trunk/libs/random/doc/Jamfile.v2	(original)
+++ trunk/libs/random/doc/Jamfile.v2	2010-06-21 00:03:55 EDT (Mon, 21 Jun 2010)
@@ -20,6 +20,7 @@
     binomial_distribution
     cauchy_distribution
     discard_block
+    discrete_distribution
     exponential_distribution
     gamma_distribution
     geometric_distribution

Modified: trunk/libs/random/doc/distributions.qbk
==============================================================================
--- trunk/libs/random/doc/distributions.qbk	(original)
+++ trunk/libs/random/doc/distributions.qbk	2010-06-21 00:03:55 EDT (Mon, 21 Jun 2010)
@@ -44,6 +44,7 @@
                              [tossing a coin 20 times and counting how many
                                 front sides are shown]]
   [[__cauchy_distribution][cauchy distribution][-]]
+  [[__discrete_distribution][discrete distribution with specific probabilities][rolling an unfair die]]
   [[__gamma_distribution][gamma distribution][-]]
   [[__poisson_distribution][poisson distribution]
                            [counting the number of alpha particles emitted

Modified: trunk/libs/random/doc/random.qbk
==============================================================================
--- trunk/libs/random/doc/random.qbk	(original)
+++ trunk/libs/random/doc/random.qbk	2010-06-21 00:03:55 EDT (Mon, 21 Jun 2010)
@@ -64,6 +64,7 @@
 [def __bernoulli_distribution [classref boost::bernoulli_distribution bernoulli_distribution]]
 [def __binomial_distribution [classref boost::random::binomial_distribution binomial_distribution]]
 [def __cauchy_distribution [classref boost::cauchy_distribution cauchy_distribution]]
+[def __discrete_distribution [classref boost::random::discrete_distribution discrete_distribution]]
 [def __gamma_distribution [classref boost::gamma_distribution gamma_distribution]]
 [def __poisson_distribution [classref boost::random::poisson_distribution poisson_distribution]]
 [def __geometric_distribution [classref boost::geometric_distribution geometric_distribution]]

Modified: trunk/libs/random/test/Jamfile.v2
==============================================================================
--- trunk/libs/random/test/Jamfile.v2	(original)
+++ trunk/libs/random/test/Jamfile.v2	2010-06-21 00:03:55 EDT (Mon, 21 Jun 2010)
@@ -45,6 +45,8 @@
 run test_binomial_distribution.cpp /boost//unit_test_framework ;
 run test_poisson.cpp ;
 run test_poisson_distribution.cpp /boost//unit_test_framework ;
+run test_discrete.cpp ;
+run test_discrete_distribution.cpp /boost//unit_test_framework ;
 
 # run nondet_random_speed.cpp ;
 # run random_device.cpp ;

Added: trunk/libs/random/test/test_discrete.cpp
==============================================================================
--- (empty file)
+++ trunk/libs/random/test/test_discrete.cpp	2010-06-21 00:03:55 EDT (Mon, 21 Jun 2010)
@@ -0,0 +1,123 @@
+/* test_poisson.cpp
+ *
+ * Copyright Steven Watanabe 2010
+ * 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)
+ *
+ * $Id$
+ *
+ */
+
+#include <boost/random/discrete_distribution.hpp>
+#include <boost/random/uniform_int.hpp>
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/lexical_cast.hpp>
+#include <boost/exception/diagnostic_information.hpp>
+#include <vector>
+#include <iostream>
+#include <numeric>
+
+#include "chi_squared_test.hpp"
+
+bool do_test(int n, long long max) {
+    std::cout << "running discrete(p0, p1, ..., p" << n-1 << ")" << " " << max << " times: " << std::flush;
+
+    std::vector<double> expected;
+    {
+        boost::mt19937 egen;
+        for(int i = 0; i < n; ++i) {
+            expected.push_back(egen());
+        }
+        double sum = std::accumulate(expected.begin(), expected.end(), 0.0);
+        for(std::vector<double>::iterator iter = expected.begin(), end = expected.end(); iter != end; ++iter) {
+            *iter /= sum;
+        }
+    }
+    
+    boost::random::discrete_distribution<> dist(expected);
+    boost::mt19937 gen;
+    std::vector<long long> results(expected.size());
+    for(long long i = 0; i < max; ++i) {
+        ++results[dist(gen)];
+    }
+
+    long long sum = std::accumulate(results.begin(), results.end(), 0ll);
+    if(sum != max) {
+        std::cout << "*** Failed: incorrect total: " << sum << " ***" << std::endl;
+        return false;
+    }
+    double chsqr = chi_squared_test(results, expected, max);
+
+    bool result = chsqr < 0.99;
+    const char* err = result? "" : "*";
+    std::cout << std::setprecision(17) << chsqr << err << std::endl;
+
+    std::cout << std::setprecision(6);
+
+    return result;
+}
+
+bool do_tests(int repeat, int max_n, long long trials) {
+    boost::mt19937 gen;
+    boost::uniform_int<> idist(1, max_n);
+    int errors = 0;
+    for(int i = 0; i < repeat; ++i) {
+        if(!do_test(idist(gen), trials)) {
+            ++errors;
+        }
+    }
+    if(errors != 0) {
+        std::cout << "*** " << errors << " errors detected ***" << std::endl;
+    }
+    return errors == 0;
+}
+
+int usage() {
+    std::cerr << "Usage: test_discrete -r <repeat> -n <max n> -t <trials>" << std::endl;
+    return 2;
+}
+
+template<class T>
+bool handle_option(int& argc, char**& argv, char opt, T& value) {
+    if(argv[0][1] == opt && argc > 1) {
+        --argc;
+        ++argv;
+        value = boost::lexical_cast<T>(argv[0]);
+        return true;
+    } else {
+        return false;
+    }
+}
+
+int main(int argc, char** argv) {
+    int repeat = 10;
+    int max_n = 100000;
+    long long trials = 1000000ll;
+
+    if(argc > 0) {
+        --argc;
+        ++argv;
+    }
+    while(argc > 0) {
+        if(argv[0][0] != '-') return usage();
+        else if(!handle_option(argc, argv, 'r', repeat)
+             && !handle_option(argc, argv, 'n', max_n)
+             && !handle_option(argc, argv, 't', trials)) {
+            return usage();
+        }
+        --argc;
+        ++argv;
+    }
+
+    try {
+        if(do_tests(repeat, max_n, trials)) {
+            return 0;
+        } else {
+            return EXIT_FAILURE;
+        }
+    } catch(...) {
+        std::cerr << boost::current_exception_diagnostic_information() << std::endl;
+        return EXIT_FAILURE;
+    }
+}

Added: trunk/libs/random/test/test_discrete_distribution.cpp
==============================================================================
--- (empty file)
+++ trunk/libs/random/test/test_discrete_distribution.cpp	2010-06-21 00:03:55 EDT (Mon, 21 Jun 2010)
@@ -0,0 +1,163 @@
+/* test_discrete_distribution.cpp
+ *
+ * Copyright Steven Watanabe 2010
+ * 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)
+ *
+ * $Id$
+ *
+ */
+
+#include <boost/random/discrete_distribution.hpp>
+#include <boost/random/linear_congruential.hpp>
+#include <boost/assign/list_of.hpp>
+#include <sstream>
+#include <vector>
+
+#define BOOST_TEST_MAIN
+#include <boost/test/unit_test.hpp>
+
+struct gen {
+    double operator()(double arg) {
+        if(arg < 100) return 100;
+        else if(arg < 103) return 1;
+        else if(arg < 107) return 2;
+        else if(arg < 111) return 1;
+        else if(arg < 114) return 4;
+        else return 100;
+    }
+};
+
+#define CHECK_PROBABILITIES(actual, expected)       \
+    do {                                            \
+        std::vector<double> _actual = (actual);     \
+        std::vector<double> _expected = (expected); \
+        BOOST_CHECK_EQUAL_COLLECTIONS(              \
+            _actual.begin(), _actual.end(),         \
+            _expected.begin(), _expected.end());    \
+    } while(false)
+
+using boost::assign::list_of;
+
+BOOST_AUTO_TEST_CASE(test_constructors) {
+    boost::random::discrete_distribution<> dist;
+    CHECK_PROBABILITIES(dist.probabilities(), list_of(1.0));
+
+#ifndef BOOST_NO_INITIALIZER_LISTS
+    boost::random::discrete_distribution<> dist_il = { 1, 2, 1, 4 };
+    CHECK_PROBABILITIES(dist_il.probabilities(), list_of(.125)(.25)(.125)(.5));
+#endif
+    std::vector<double> probs = boost::assign::list_of(1.0)(2.0)(1.0)(4.0);
+
+    boost::random::discrete_distribution<> dist_r(probs);
+    CHECK_PROBABILITIES(dist_r.probabilities(), list_of(.125)(.25)(.125)(.5));
+    
+    boost::random::discrete_distribution<> dist_it(probs.begin(), probs.end());
+    CHECK_PROBABILITIES(dist_it.probabilities(), list_of(.125)(.25)(.125)(.5));
+    
+    boost::random::discrete_distribution<> dist_fun(4, 99, 115, gen());
+    CHECK_PROBABILITIES(dist_fun.probabilities(), list_of(.125)(.25)(.125)(.5));
+
+    boost::random::discrete_distribution<> copy(dist);
+    BOOST_CHECK_EQUAL(dist, copy);
+    boost::random::discrete_distribution<> copy_r(dist_r);
+    BOOST_CHECK_EQUAL(dist_r, copy_r);
+
+    boost::random::discrete_distribution<> notpow2(3, 99, 111, gen());
+    BOOST_REQUIRE_EQUAL(notpow2.probabilities().size(), 3u);
+    BOOST_CHECK_CLOSE_FRACTION(notpow2.probabilities()[0], 0.25, 0.00000000001);
+    BOOST_CHECK_CLOSE_FRACTION(notpow2.probabilities()[1], 0.50, 0.00000000001);
+    BOOST_CHECK_CLOSE_FRACTION(notpow2.probabilities()[2], 0.25, 0.00000000001);
+    boost::random::discrete_distribution<> copy_notpow2(notpow2);
+    BOOST_CHECK_EQUAL(notpow2, copy_notpow2);
+}
+
+BOOST_AUTO_TEST_CASE(test_param) {
+    std::vector<double> probs = boost::assign::list_of(1.0)(2.0)(1.0)(4.0);
+    boost::random::discrete_distribution<> dist(probs);
+    boost::random::discrete_distribution<>::param_type param = dist.param();
+    CHECK_PROBABILITIES(param.probabilities(), list_of(.125)(.25)(.125)(.5));
+    boost::random::discrete_distribution<> copy1(param);
+    BOOST_CHECK_EQUAL(dist, copy1);
+    boost::random::discrete_distribution<> copy2;
+    copy2.param(param);
+    BOOST_CHECK_EQUAL(dist, copy2);
+
+    boost::random::discrete_distribution<>::param_type param_copy = param;
+    BOOST_CHECK_EQUAL(param, param_copy);
+    BOOST_CHECK(param == param_copy);
+    BOOST_CHECK(!(param != param_copy));
+    boost::random::discrete_distribution<>::param_type param_default;
+    CHECK_PROBABILITIES(param_default.probabilities(), list_of(1.0));
+    BOOST_CHECK(param != param_default);
+    BOOST_CHECK(!(param == param_default));
+    
+#ifndef BOOST_NO_INITIALIZER_LISTS
+    boost::random::discrete_distribution<>::param_type
+        parm_il = { 1, 2, 1, 4 };
+    CHECK_PROBABILITIES(parm_il.probabilities(), list_of(.125)(.25)(.125)(.5));
+#endif
+
+    boost::random::discrete_distribution<>::param_type parm_r(probs);
+    CHECK_PROBABILITIES(parm_r.probabilities(), list_of(.125)(.25)(.125)(.5));
+    
+    boost::random::discrete_distribution<>::param_type
+        parm_it(probs.begin(), probs.end());
+    CHECK_PROBABILITIES(parm_it.probabilities(), list_of(.125)(.25)(.125)(.5));
+    
+    boost::random::discrete_distribution<>::param_type
+        parm_fun(4, 99, 115, gen());
+    CHECK_PROBABILITIES(parm_fun.probabilities(), list_of(.125)(.25)(.125)(.5));
+}
+
+BOOST_AUTO_TEST_CASE(test_min_max) {
+    std::vector<double> probs = boost::assign::list_of(1.0)(2.0)(1.0);
+    boost::random::discrete_distribution<> dist;
+    BOOST_CHECK_EQUAL((dist.min)(), 0);
+    BOOST_CHECK_EQUAL((dist.max)(), 0);
+    boost::random::discrete_distribution<> dist_r(probs);
+    BOOST_CHECK_EQUAL((dist_r.min)(), 0);
+    BOOST_CHECK_EQUAL((dist_r.max)(), 2);
+}
+
+BOOST_AUTO_TEST_CASE(test_comparison) {
+    std::vector<double> probs = boost::assign::list_of(1.0)(2.0)(1.0)(4.0);
+    boost::random::discrete_distribution<> dist;
+    boost::random::discrete_distribution<> dist_copy(dist);
+    boost::random::discrete_distribution<> dist_r(probs);
+    boost::random::discrete_distribution<> dist_r_copy(dist_r);
+    BOOST_CHECK(dist == dist_copy);
+    BOOST_CHECK(!(dist != dist_copy));
+    BOOST_CHECK(dist_r == dist_r_copy);
+    BOOST_CHECK(!(dist_r != dist_r_copy));
+    BOOST_CHECK(dist != dist_r);
+    BOOST_CHECK(!(dist == dist_r));
+}
+
+BOOST_AUTO_TEST_CASE(test_streaming) {
+    std::vector<double> probs = boost::assign::list_of(1.0)(2.0)(1.0)(4.0);
+    boost::random::discrete_distribution<> dist(probs);
+    std::stringstream stream;
+    stream << dist;
+    boost::random::discrete_distribution<> restored_dist;
+    stream >> restored_dist;
+    BOOST_CHECK_EQUAL(dist, restored_dist);
+}
+
+BOOST_AUTO_TEST_CASE(test_generation) {
+    std::vector<double> probs = boost::assign::list_of(0.0)(1.0);
+    boost::minstd_rand0 gen;
+    boost::random::discrete_distribution<> dist;
+    boost::random::discrete_distribution<> dist_r(probs);
+    for(int i = 0; i < 10; ++i) {
+        int value = dist(gen);
+        BOOST_CHECK_EQUAL(value, 0);
+        int value_r = dist_r(gen);
+        BOOST_CHECK_EQUAL(value_r, 1);
+        int value_param = dist_r(gen, dist.param());
+        BOOST_CHECK_EQUAL(value_param, 0);
+        int value_r_param = dist(gen, dist_r.param());
+        BOOST_CHECK_EQUAL(value_r_param, 1);
+    }
+}