From 1c4bce4f96178fb8806ff1eec030253df4ce8fe1 Mon Sep 17 00:00:00 2001 From: bookworm6 Date: Thu, 2 Jul 2026 19:58:51 -0700 Subject: [PATCH] Added dynamic_discrete_distribution dynamic_discrete_distribution is a dynamic_distribution with updatable weights. Selection and update have been highly optimized --- doc/Jamfile.v2 | 1 + doc/distributions.qbk | 1 + doc/random.qbk | 1 + include/boost/random.hpp | 1 + .../random/dynamic_discrete_distribution.hpp | 849 ++++++++++++++++++ test/Jamfile.v2 | 1 + test/multiprecision_int_test.cpp | 5 + test/test_dynamic_discrete_distribution.cpp | 347 +++++++ 8 files changed, 1206 insertions(+) create mode 100644 include/boost/random/dynamic_discrete_distribution.hpp create mode 100644 test/test_dynamic_discrete_distribution.cpp diff --git a/doc/Jamfile.v2 b/doc/Jamfile.v2 index 7a7bc91c5..8d069ab22 100644 --- a/doc/Jamfile.v2 +++ b/doc/Jamfile.v2 @@ -25,6 +25,7 @@ doxygen_files = chi_squared_distribution discard_block discrete_distribution + dynamic_discrete_distribution exponential_distribution extreme_value_distribution faure diff --git a/doc/distributions.qbk b/doc/distributions.qbk index 7c5a4f948..94b1c5763 100644 --- a/doc/distributions.qbk +++ b/doc/distributions.qbk @@ -95,6 +95,7 @@ statistically to it are not acceptable. [table Sampling Distributions [[distribution] [explanation] [example]] [[__discrete_distribution][discrete distribution with specific probabilities][rolling an unfair die]] + [[__dynamic_discrete_distribution][discrete distribution with specific updatable probabilities][rolling an unfair die and changing the weights]] [[__piecewise_constant_distribution][-][-]] [[__piecewise_linear_distribution][-][-]] ] diff --git a/doc/random.qbk b/doc/random.qbk index cbc5f91ec..0e2a451fd 100644 --- a/doc/random.qbk +++ b/doc/random.qbk @@ -89,6 +89,7 @@ [def __binomial_distribution [classref boost::random::binomial_distribution binomial_distribution]] [def __cauchy_distribution [classref boost::random::cauchy_distribution cauchy_distribution]] [def __discrete_distribution [classref boost::random::discrete_distribution discrete_distribution]] +[def __dynamic_discrete_distribution [classref boost::random::dynamic_discrete_distribution dynamic_discrete_distribution]] [def __gamma_distribution [classref boost::random::gamma_distribution gamma_distribution]] [def __hyperexponential_distribution [classref boost::random::hyperexponential_distribution hyperexponential_distribution]] [def __laplace_distribution [classref boost::random::laplace_distribution laplace_distribution]] diff --git a/include/boost/random.hpp b/include/boost/random.hpp index 884aa93ec..bcd149beb 100644 --- a/include/boost/random.hpp +++ b/include/boost/random.hpp @@ -64,6 +64,7 @@ #include #include #include +#include #include #include #include diff --git a/include/boost/random/dynamic_discrete_distribution.hpp b/include/boost/random/dynamic_discrete_distribution.hpp new file mode 100644 index 000000000..5b8f35647 --- /dev/null +++ b/include/boost/random/dynamic_discrete_distribution.hpp @@ -0,0 +1,849 @@ +#ifndef BOOST_RANDOM_DYNAMIC_DISCRETE_DISTRIBUTION_HPP_INCLUDED +#define BOOST_RANDOM_DYNAMIC_DISCRETE_DISTRIBUTION_HPP_INCLUDED +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + +namespace boost { namespace random{ + +///////////////////////// +// Underlying tree data structure +namespace detail{ +template + class complete_kary_complete_tree { + static constexpr typename std::make_unsigned::type unsigned_fanout = Fanout; + + + public: + + // Compute minimal complete k-ary tree size to store exactly n leaves + // without needing bounds checks during selection. returns the index of the first leaf and the total number of leaves + // Returns {total_nodes_excluding_root, leaf_start_index} + std::pair minimal_tree_shape(LiteralIntType n) { + + if (n == 0) return {0, 0}; // no internal nodes, no leaves + + LiteralIntType k = boost::core::countr_zero(unsigned_fanout); + typename std::make_unsigned::type unsigned_n = static_cast::type>(n); + LiteralIntType max_leaf = LiteralIntType(1) << (((boost::core::bit_width(unsigned_n - 1) + k - 1) / k) * k); + + LiteralIntType leaf_start = 0; //the index of the first node that can contain leaves + LiteralIntType level = 0; //the index of the first node in the bottom row + + while (true) { + level = first_child_of(level); + if (level >= max_leaf) break; + leaf_start=level; + } + + LiteralIntType total_nodes = leaf_start+n; + return {total_nodes, leaf_start}; + } + + + + static_assert((Fanout & (Fanout - 1)) == 0, "fanout must be power of two"); + + using position_type = LiteralIntType; + + //constructors + complete_kary_complete_tree() : data_(1, Real(0)), leaf_start_(0), leaf_count_(0) {} + + explicit complete_kary_complete_tree(LiteralIntType leaf_count) { + resize(leaf_count); + } + + // Resize to accommodate exactly n leaves + void resize(LiteralIntType n) { + if (n == 0) { + data_.assign(1, Real(0)); + leaf_start_ = 0; + leaf_count_ = 0; + return; + } + + std::pair treeshape = minimal_tree_shape(n); + LiteralIntType total_nodes = treeshape.first; + // allocate all internal nodes + leaves (root is kept separately by the derived class) + data_.resize(total_nodes); + leaf_count_ = n; + leaf_start_ = treeshape.second; + } + + + + LiteralIntType size() const { + return data_.size(); + } + + Real& value_of(position_type p) { + BOOST_ASSERT(p < data_.size()); + + return data_[p]; + } + + const Real& value_of(position_type p) const { + BOOST_ASSERT(p < data_.size()); + return data_[p]; + } + + // Tree navigation (-1 based indexing). root stored separately + + + position_type parent_of(position_type i) const { + BOOST_ASSERT(i > 0); + return (i >> log2_fanout) - 1; + } + + position_type first_child_of(position_type i) const { + return (i + 1) << log2_fanout; + } + + position_type child_index(position_type parent, LiteralIntType child_num) const { + return first_child_of(parent) + child_num; + } + + bool is_leaf(position_type i) const { + return i >= leaf_start_; + } + + LiteralIntType leaf_start() const { + return leaf_start_; + } + + LiteralIntType leaf_count() const { + return leaf_count_; + } + + std::vector > &data() {return data_;} + + const std::vector > &data() const {return data_;} + + private: + std::vector > data_; + LiteralIntType leaf_start_; + LiteralIntType leaf_count_; + LiteralIntType max_leaf_; + static constexpr LiteralIntType fanout = Fanout; + + + + static constexpr LiteralIntType log2_fanout = core::countr_zero(unsigned_fanout); + };} + +/////////////////////// + + +/** + * @brief dynamic_discrete_distribution is a highly optimized library for weighted random selection modeled after std::discrete_distribution. It selects an index where the index's probability of being selected is weighted by inputted weights. It satisfies all RandomNumberDistribution requirements except for the requirement of constant time equality comparisons between both distribution objects and parameter objects, which is not possible in discrete distributions. + * * This library is designed to be used in discrete event simulation applications in which there is a set of events with different probabilities of occurring, and the simulation must choose an event to occur and update probabilities of other events accordingly. These simulations often involve many more updates than selections, but exact ratios of update to selection differ depending on the application. This library provides efficient update and selection in which the underlying tree data structure can be tuned at compile time to prioritize update over selection to varying degrees. + * @tparam IntType is the type of the integers returned by operator(), which is the selection function. IntType must be an unsigned integer type + * @tparam Real is the type of the weights + * @tparam Fanout, which must be an positive integer power of 2, controls the branching factor of the underlying complete tree data structure and can be adjusted to change the amount that update is prioritized over selection. update_weight has a runtime of $O(log_{\text{fanout}} N)$ while selection (operator ()) has a runtime of $O( \text{fanout} * (log_{\text{fanout}} N))$ . Additionally, trees with larger fanouts use less memory. + * @tparam Precision controls the number of bits of randomness generated during selection. + */ +template < + class IntType = int, + class Real = double, + int Fanout = 16, + size_t Precision = std::numeric_limits::digits +> +class dynamic_discrete_distribution { + //This protects against non literal IntType types (necessary because of the boost::multiprecision::cpp_int types). indexes can't be larger than size_t anyway, so we just use size_t. We use this in the place of IntType for internal work, and cast to IntType before returning + using LiteralIntType = typename std::conditional::value && std::is_constructible::value,IntType,size_t>::type; + + static constexpr typename std::make_unsigned::type unsigned_fanout = Fanout; + static_assert(Fanout>0 && boost::core::has_single_bit(unsigned_fanout),"template parameter Fanout must be a positive power of 2"); + using This = dynamic_discrete_distribution; + + +public: + using input_type = Real; + using result_type = IntType; + + /** + * @brief standard library random number distributions separate state related to the parameters of a distribution and state related to the generation of random numbers by defining a member class param_type that holds the distribution parameters and related data structures. The distribution stores its parameter set in a param_type object. + */ + class Param : protected detail::complete_kary_complete_tree { + using base_tree = detail::complete_kary_complete_tree; + using pos_type = typename base_tree::position_type; + + public: + using distribution_type = This; + + /** + * @brief constructs a `param_type` parameter set with a single element 0 with a weight of 1. + */ + Param() : Param({1.0}) {} + + /** + * @brief constructs a `param_type` parameter set with weights equal to the value in `weights` at each index. + */ + Param(const std::vector& weights) + : Param(weights.begin(), weights.end()) {} + + /** + * @brief constructs a `param_type` parameter set with the weights in the initializer_list `il`. Each weight will have the same index as it did in il. + */ + Param(const std::initializer_list& il) + : Param(il.begin(), il.end()) {} + + /** + * @brief Constructs a `param_type` parameter set with `count` weights that are generated using function `unary_op`. Each of the weights is equal to ${w_i} = unary_op(xmin + δ(i + 0.5))$, where $δ = (xmax − xmin)count$ and $i ∈ {0, ..., count − 1}$. `xmin` and `xmax` must be such that `δ > 0`. If `count == 0` the effects are the same as of the default constructor. + */ + template< class UnaryOperation > + Param( std::size_t count, double xmin, double xmax, UnaryOperation unary_op ) + : Param(weight_distribution(count,xmin,xmax,unary_op)) {} + + /** + * @brief constructs the distribution from the first and last iterators to a collection of weights of type `weight_type`. Each weight is indexed according to the number of elements between itself and first. + */ + template + Param(InputIt first, InputIt last) + : base_tree() + { + size_t distance = std::distance(first, last); + LiteralIntType n = distance; + typename std::make_unsigned::type n_unsigned = n; + + //n = 2; + + // Round up leaves to nearest full complete k-ary tree level (power of fanout) + if (n <= fanout){ + max_leaf_ = fanout; + num_layers_ = 1; + } + else{ + LiteralIntType k = boost::core::countr_zero(unsigned_fanout); + num_layers_ = (boost::core::bit_width(n_unsigned - 1) + k - 1) / k; + max_leaf_ = LiteralIntType(1) << (((boost::core::bit_width(n_unsigned - 1) + k - 1) / k) * k); + } + + leaf_start_ = base_tree::minimal_tree_shape(n).second; + base_tree::resize(n); + leaf_end_ = n; + + + // Copy weights to leaves, pad with zeros + InputIt it = first; + for (LiteralIntType i = leaf_start_; i < leaf_start_ + n; ++i) { + weightsum_of(i) = std::max(Real(*it), Real(0)); + ++it; + } + + // Build sums bottom-up from leaves to root (excluding root) + for (ptrdiff_t i = leaf_start_ - 1; i >= 0; --i) { + Real sum = 0; + pos_type first_child = (i + 1) * fanout; // -1-index adjustment + for (LiteralIntType c = 0; c < fanout; ++c) { + pos_type child = first_child + c; + if (child >= leaf_start_+n) break; + sum += weightsum_of(child); + } + weightsum_of(i) = sum; + } + + // Compute root separately + Real sum = 0; + pos_type first_child = 0 * fanout; // root's first child in array + if (leaf_start_ == 0){ + for (LiteralIntType c = 0; c < leaf_end_; ++c) { + pos_type child = first_child + c; + if (child >= max_leaf_) break; + sum += weightsum_of(child); + } + } + else{ + for (LiteralIntType c = 0; c < fanout; ++c) { + pos_type child = first_child + c; + if (child >= max_leaf_) break; + sum += weightsum_of(child); + } + } + total_weight_ = sum; + } + + /** + * @brief returns a vector of probabilities of each integer that could be generated by a dynamic_discrete_distribution using this parameter set. + */ + std::vector probabilities() const { + std::vector probs(leaf_end_, Real(0)); + Real total = total_weight(); + if (total <= Real(0)) return probs; + + for (LiteralIntType i = 0; i < leaf_end_; ++i) { + probs[i] = weightsum_of(leaf_start_ + i) / total; + } + return probs; + } + + /** + * @brief returns the minimum integer that could be generated by a dynamic_discrete_distribution using this parameter set assuming that it is not empty. This minimum will always be 0. + */ + static constexpr result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () { return 0; } + + + /** + * @brief returns the maximum integer that could be generated by a dynamic_discrete_distribution using this parameter set assuming that it is not empty. + */ + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return leaf_end_ == 0 ? 0 : leaf_end_ - 1; } + + /** + * @brief updates weight of int`i` to `new_weight`. If i is not in this parameter set, the behavior is undefined. + */ + void update_weight(IntType i, Real new_weight) { + LiteralIntType i_literal = i; + BOOST_ASSERT(new_weight >= Real(0)); + BOOST_ASSERT(i_literal <= leaf_end_); + BOOST_ASSERT(i_literal>=0); + i_literal = leaf_start_ + i_literal; + Real diff = new_weight - weightsum_of(i_literal); + weightsum_of(i_literal) = new_weight; + total_weight_ += diff; + + + while (i_literal >=fanout) { + i_literal = base_tree::parent_of(i_literal); + + weightsum_of(i_literal) += diff; + } + } + + /** + * @brief gets the weight of int `i`. If i is not in the parameter set, the behavior is undefined. + */ + Real get_weight(IntType i) const { + return get_weight_literal(i); + } + + /** + * @brief returns the number of elements in this parameter set. + */ + IntType size() const { return leaf_end_; } + + /** + * @brief returns the total of all of the weights in this parameter set. + */ + Real total_weight() const noexcept { + if (base_tree::size() == 0) return Real(0); + return total_weight_; + } + + /** + * @brief adds weights.size() elements to the end of this parameter set with weights given by `weights` vector. + */ + void push_back(const std::vector& weights){ + for (Real w : weights) { + push_back(w); + } + } + + /** + * @brief adds an element to the end with weight `weight`. + */ + void push_back(Real weight) { + expand(leaf_end_+1); + update_weight(leaf_end_, weight); + leaf_end_++; + } + + /** + * @brief removes the `count` highest integers. If the parameter set is empty, this function's behavior is undefined. NOTE: in this implementation, this does not deallocate memory or decrease the depth of the underlying tree based data structure. For the purposes of asymptotic analysis, consider "N" to be the largest number of elements ever held in this parameter set. + */ + void pop_back(IntType count) { + for(IntType i = 0; i < count; ++i) { + pop_back(); + } + } + + /** + * @brief removes the highest integer. If the parameter set is empty, this function's behavior is undefined. NOTE: in this implementation, this does not deallocate memory or decrease the depth of the underlying tree based data structure. For the purposes of asymptotic analysis, consider "N" to be the largest number of elements ever held in this parameter set. + */ + void pop_back() { + BOOST_ASSERT(leaf_end_ > 0); + leaf_end_--; + update_weight(leaf_end_, Real(0)); + } + + /** + * @brief compares the weights of each element in `rhs` and `lhs` for equality. + */ + BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(Param, lhs, rhs){ + LiteralIntType sizer = rhs.literal_size(); + if (sizer != lhs.literal_size()){ + return false; + } + for (LiteralIntType i=0; i q; // use signed for -1 root +// q.push(-1); // -1 represents the root stored separately + +// IntType level = 0; + +// while (!q.empty()) { +// IntType level_size = q.size(); +// os << "Level " << level << ": "; + +// for (IntType i = 0; i < level_size; ++i) { +// ptrdiff_t node = q.front(); +// q.pop(); + +// Real w = (node == -1) ? total_weight_ : weightsum_of(static_cast(node)); +// os << "[" << node << "]=" << w << " "; + +// // enqueue children +// pos_type first_child = (node + 1) * fanout; // -1-index adjustment +// for (IntType c = 0; c < fanout; ++c) { +// pos_type child = first_child + c; +// if (child >= total_nodes) break; +// q.push(child); +// } +// } + +// os << "\n"; +// ++level; +// } +// } + + /** + * @brief compares the weights of each element in `rhs` and `lhs` for inequality + */ + BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(Param) + + /** + * @brief Restores the parameters with data read from `stream`. The formatting flags of `stream` are unchanged. The data must have been written using a stream with the same locale, `CharT` and `Traits` template parameters, otherwise the behavior is undefined. If bad input is encountered, stream.setstate(std::ios::failbit) is called, which may throw std::ios_base::failure. the parameter set is unchanged in that case. + */ + BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(stream, Param, dist){ + std::vector newWeights; + detail::read_vector(stream, newWeights); + if (stream){ + dist = Param(newWeights); + } + + + return stream; + } + + /** + * @brief Writes a textual representation of the parameters to `stream`. The formatting flags and fill character of `stream` are unchanged. + */ + BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(stream, Param, dist){ + std::vector leaves(dist.base_tree::data().begin()+dist.leaf_start_, dist.base_tree::data().begin()+dist.leaf_start_+dist.leaf_end_ ); + detail::print_vector(stream, leaves); + + return stream; + } + + private: + //size for internal use, involving fewer casts. + LiteralIntType literal_size()const{ + return leaf_end_; + } + + //get_weight for internal use + Real get_weight_literal(LiteralIntType i) const { + BOOST_ASSERT(i <= leaf_end_+leaf_start_); + BOOST_ASSERT(i>=0); + return weightsum_of(leaf_start_ + i); + } + + Real& weightsum_of(pos_type p) { return base_tree::value_of(p); } + const Real& weightsum_of(pos_type p) const { return base_tree::value_of(p); } + + // Expand from current leaf_count_ to new_leaf_count (must be larger), structurally, without recomputing + // new expand: argument is NEW_LEAF_COUNT (number of leaves you want after expansion) + void expand(LiteralIntType new_leaf_count) { + if (new_leaf_count <= leaf_end_) return; // nothing to do + auto &data_ = base_tree::data(); + if (new_leaf_count > max_leaf_) { + max_leaf_ *= fanout; + // old shape + std::pair oldShape = base_tree::minimal_tree_shape(leaf_end_); + LiteralIntType old_total_nodes = oldShape.first; + + // new shape + std::pair newShape = base_tree::minimal_tree_shape(new_leaf_count); + LiteralIntType new_total_nodes = newShape.first; + leaf_start_ = newShape.second; + + std::vector old_starts; + std::vector old_levels; + + std::vector new_starts; + std::vector new_levels; + + LiteralIntType position = 0; + LiteralIntType sum = 0; + + // Creating the level lengths and level starts for the old size + while (true) { + position = base_tree::first_child_of(position); + old_starts.push_back(sum); + LiteralIntType level = position - sum; + sum += level; + old_levels.push_back(level); + if (position >= old_total_nodes) break; + } + + // Creating the level lengths and level starts for the new size + position = 0; + sum = 0; + while (true) { + position = base_tree::first_child_of(position); + new_starts.push_back(sum); + LiteralIntType level = position - sum; + sum += level; + new_levels.push_back(level); + if (position >= new_total_nodes) break; + } + + // allocate new array (zero-initialized) + std::vector> new_data(new_total_nodes, Real(0)); + + // copy each old level block into the next-deeper level of the new layout. + // old level i -> new level (i+1). That packs the old block contiguously at the start + // of the larger new level (the remainder stays zero). + for (LiteralIntType i = 0; i < old_levels.size(); ++i) { + LiteralIntType old_start = old_starts[i]; + LiteralIntType old_sz = old_levels[i]; + LiteralIntType new_level_index = i + 1; // destination level index + LiteralIntType new_start = new_starts[new_level_index]; + + + std::copy_n(data_.begin() + old_start, old_sz, new_data.begin() + new_start); + } + + // recompute only the new top internal layer (level 0) from its children (level 1) + // top-level nodes occupy global indices new_starts[0] .. new_starts[0]+new_levels[0]-1 + // their first child indices can be computed with first_child_of(parent_index) + LiteralIntType top_count = new_levels[0]; + LiteralIntType top_start = new_starts[0]; // usually 0 + for (LiteralIntType j = 0; j < top_count; ++j) { + LiteralIntType parent_idx = top_start + j; + // first child in global indexing: + LiteralIntType first_child = base_tree::first_child_of(parent_idx); + Real sum = Real(0); + for (LiteralIntType c = 0; c < fanout; ++c) { + LiteralIntType child = first_child + c; + if (child >= new_total_nodes) break; + sum += new_data[child]; + } + new_data[parent_idx] = sum; + } + // commit! + num_layers_++; + data_.swap(new_data); + + } + else if (new_leaf_count+leaf_start_ >= base_tree::size()) { + base_tree::resize(std::min(new_leaf_count*2,leaf_start_+max_leaf_)); + } + } + + LiteralIntType max_leaf_; + LiteralIntType leaf_end_; // number of leaves requested by user + typename std::make_signed::type leaf_start_; // index of first leaf in data_ + LiteralIntType num_layers_; + Real total_weight_ = 0; + + template< class UnaryOperation > + static std::vector weight_distribution(std::size_t count, double xmin, double xmax,UnaryOperation unary_op ){ + std::vector distro; + if (count<=0){ + distro.push_back(1); + } + else{ + double delta = (xmax - xmin)/count; + for (IntType i=0; i& weights) + : dynamic_discrete_distribution(weights.begin(), weights.end()) {} + + /** + * @brief constructs a distribution with the parameters in p. + */ + explicit dynamic_discrete_distribution(const Param& p){ + param_ = p; + } + + /** + * @brief constructs a distribution with the weights in the initializer_list `il`. Each weight will have the same index as it did in il. + */ + dynamic_discrete_distribution(const std::initializer_list& il) + : dynamic_discrete_distribution(il.begin(), il.end()) {} + + /** + * @brief Constructs the distribution with `count` weights that are generated using function `unary_op`. Each of the weights is equal to ${w_i} = unary_op(xmin + δ(i + 0.5))$, where $δ = (xmax − xmin)count$ and $i ∈ {0, ..., count − 1}$. `xmin` and `xmax` must be such that `δ > 0`. If `count == 0` the effects are the same as of the default constructor. + */ + template< class UnaryOperation > + dynamic_discrete_distribution( std::size_t count, double xmin, double xmax, UnaryOperation unary_op ){ + param_ = Param(count,xmin,xmax,unary_op); + } + + /** + * @brief constructs the distribution from the first and last iterators to a collection of weights of type `weight_type`. Each weight is indexed according to the number of elements between itself and first. + */ + template + dynamic_discrete_distribution(InputIt first, InputIt last){ + param_ = Param(first,last); + } + + /** + * @brief Resets the internal state of the distribution object. After a call to this function, the next call to operator() on the distribution object will not be dependent on previous calls to operator(). The distribution still depends on past calls to update_weight, push, and pop. (Note: in this implementation there is no need for reset, so it does no work. It is present for compatibility reasons) + */ + void reset() {} + + + /** + * @brief Generates random numbers distributed according to the weights in the parameter set `param`. If it does not contain at least 1 weight, the behavior is undefined + */ + template + IntType operator()(URNG& g,const Param& param) const { + Real total = param.total_weight(); + assert(total>0 && param.leaf_end_>0); + if (total <= Real(0)) return 0; + Real target = boost::random::generate_canonical(g) * total; + if (target == Real(0)) return 0; + typename Param::pos_type first_child = 0; + + // Start at the top internal node (index 0) + typename Param::pos_type node = 0; + + for(LiteralIntType i=0; i= param.leaf_start_+param.leaf_end_){ //floating point errors made selection move into empty part of tree at some point + return this->operator()(g,param); //run selection algorithm again (probability of floating point errors causing this is quite small) + } + Real cumulative = 0; + bool chosen = false; + for (LiteralIntType c = 0; (c < fanout)&&(c+first_child= param.leaf_end_+param.leaf_start_) break; + Real w = param.weightsum_of(child); + if (target < cumulative + w) { + node = child; + target -=cumulative; + chosen = true; + + break; + } + cumulative += w; + } + if (chosen == false){ + node = first_child; + target -= cumulative; + target +=param.weightsum_of(first_child); + } + + // otherwise, node has been updated to chosen_child + + + // node is now a leaf + return static_cast(node - param.leaf_start_); + + } + /** + * @brief Generates random numbers that are distributed according to the weights in the distribution's parameter set. If it does not contain at least 1 weight, the behavior is undefined + */ + template + result_type operator()(URNG& g) const { + return this->operator()(g,param_); + } + + /** + * @brief returns a vector of probabilities of each integer that could be generated by this distribution using its associated parameter set. + */ + std::vector probabilities() const { + return param_.probabilities(); + } + + /** + * @brief gets the distribution's parameter set + */ + Param param() const { + return param_; + } + + /** + * @brief sets the distribution's parameter set to `p` + */ + void param(const Param& p) { + param_ = p; + } + + /** + * @brief returns the minimum integer that could be generated by the operator called on the distribution's parameter set assuming that it is not empty. This minimum will always be 0. + */ + static constexpr result_type min() { return 0; } + + /** + * @brief returns the maximum integer that could be generated by the operator called on the distribution's parameter set assuming that it is not empty. + */ + result_type max() const { return param_.max(); } + + /** + * @brief updates weight of int`i` in the distribution's parameter set to `new_weight`. If i is not in the distribution's parameter set, the behavior is undefined. + */ + void update_weight(IntType i, Real new_weight) { + param_.update_weight(i, new_weight); + } + + /** + * @brief gets the weight of int `i` in the distribution's parameter set. If i is not in the parameter set, the behavior is undefined. + */ + Real get_weight(IntType i) const { + return param_.get_weight_literal(i); + } + + /** + * @brief returns the number of weights in the distribution's parameter set. + */ + IntType size() const { return param_.size(); } + + /** + * @brief returns the total of all of the weights in the distribution's parameter set. + */ + Real total_weight() const noexcept { + return param_.total_weight(); + } + + /** + * @brief adds weights.size() elements to the end of the distribution's parameter set with weights given by `weights` vector. + */ + void push_back(const std::vector& weights){ + param_.push_back(weights); + } + + /** + * @brief adds an element to the end of distribution's parameter set with weight `weight`. + */ + void push_back(Real weight) { + param_.push_back(weight); + } + + /** + * @brief removes the `count` highest integers from the distribution's parameter set. If the distribution's parameter set is empty, this function's behavior is undefined. NOTE: in this implementation, this does not deallocate memory or decrease the depth of the underlying tree based data structure. For the purposes of asymptotic analysis, consider "N" to be the largest number of elements ever held in the distribution's parameter set. + */ + void pop_back(IntType count) { + param_.pop_back(count); + } + + /** + * @brief removes the highest integer from the distribution's parameter set. If the distribution's parameter set is empty, this function's behavior is undefined. NOTE: in this implementation, this does not deallocate memory or decrease the depth of the underlying tree based data structure. For the purposes of asymptotic analysis, consider "N" to be the largest number of elements ever held in the distribution's parameter set. + */ + void pop_back() { + param_.pop_back(); + } + + /** + * @brief compares the parameter sets of `rhs` and `lhs` for equality + */ + BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(This, lhs, rhs){ + return (rhs.param_ == lhs.param_); + } + + /** + * @brief compares the parameter sets of `rhs` and `lhs` for inequality + */ + BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(This) + + + /** + * @brief Restores the distribution parameters with data read from `stream`. The formatting flags of `stream` are unchanged. The data must have been written using a stream with the same locale, `CharT` and `Traits` template parameters, otherwise the behavior is undefined. If bad input is encountered, stream.setstate(std::ios::failbit) is called, which may throw std::ios_base::failure. `dist` is unchanged in that case. + */ + BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(stream, This, dist){ + return (stream>>dist.param_); + } + + /** + * @brief Writes a textual representation of the distribution parameters to `stream`. The formatting flags and fill character of `stream` are unchanged. + */ + BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(stream, This, dist){ + return (stream< + +#endif \ No newline at end of file diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 17157527d..a6bc24659 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -118,6 +118,7 @@ run test_poisson.cpp math_test ; run test_poisson_distribution.cpp /boost/test//boost_unit_test_framework ; run test_discrete.cpp math_test ; run test_discrete_distribution.cpp /boost/assign//boost_assign /boost/test//boost_unit_test_framework ; +run test_dynamic_discrete_distribution.cpp /boost/test//boost_unit_test_framework ; run test_gamma.cpp math_test ; run test_gamma_distribution.cpp /boost/test//boost_unit_test_framework ; run test_weibull.cpp math_test ; diff --git a/test/multiprecision_int_test.cpp b/test/multiprecision_int_test.cpp index a861cca39..421b958be 100644 --- a/test/multiprecision_int_test.cpp +++ b/test/multiprecision_int_test.cpp @@ -192,6 +192,11 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(distributions, distribution_type, uniform_distribu typedef boost::mpl::list < boost::random::discrete_distribution < boost::multiprecision::cpp_int, double>, boost::random::discrete_distribution + + boost::random::dynamic_discrete_distribution + boost::random::dynamic_discrete_distribution < boost::multiprecision::cpp_int, double,16>, + + > other_distributions; diff --git a/test/test_dynamic_discrete_distribution.cpp b/test/test_dynamic_discrete_distribution.cpp new file mode 100644 index 000000000..23b7e0b95 --- /dev/null +++ b/test/test_dynamic_discrete_distribution.cpp @@ -0,0 +1,347 @@ +/* test_dynamic_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 +#include +#include +#include +#include +#include +#include +#include "concepts.hpp" +#include "chi_squared_test.hpp" + +#define BOOST_TEST_MAIN +#include + + + +using boost::random::test::RandomNumberDistribution; +using boost::random::dynamic_discrete_distribution; + +BOOST_CONCEPT_ASSERT((RandomNumberDistribution< dynamic_discrete_distribution<> >)); + +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 _actual = (actual); \ + std::vector _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::dynamic_discrete_distribution<> dist; + CHECK_PROBABILITIES(dist.probabilities(),(std::vector{1.0})); + +#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST + boost::random::dynamic_discrete_distribution<> dist_il = { 1, 2, 1, 4 }; + CHECK_PROBABILITIES(dist_il.probabilities(), (std::vector{.125,.25,.125,.5})); +#endif + std::vector probs = {1.0,2.0,1.0,4.0}; + + boost::random::dynamic_discrete_distribution<> dist_r(probs); + CHECK_PROBABILITIES(dist_r.probabilities(), (std::vector{.125,.25,.125,.5})); + + boost::random::dynamic_discrete_distribution<> dist_i({1.0,2.0,1.0,4.0}); + CHECK_PROBABILITIES(dist_r.probabilities(), (std::vector{.125,.25,.125,.5})); + + boost::random::dynamic_discrete_distribution<> dist_it(probs.begin(), probs.end()); + CHECK_PROBABILITIES(dist_it.probabilities(), (std::vector{.125,.25,.125,.5})); + + boost::random::dynamic_discrete_distribution<> dist_fun(4, 99, 115, gen()); + CHECK_PROBABILITIES(dist_fun.probabilities(), (std::vector{.125,.25,.125,.5})); + boost::random::dynamic_discrete_distribution<> copy(dist); + BOOST_CHECK_EQUAL(dist, copy); + boost::random::dynamic_discrete_distribution<> copy_r(dist_r); + BOOST_CHECK_EQUAL(dist_r, copy_r); + + boost::random::dynamic_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::dynamic_discrete_distribution<> copy_notpow2(notpow2); + BOOST_CHECK_EQUAL(notpow2, copy_notpow2); +} + +BOOST_AUTO_TEST_CASE(test_param) { + std::vector probs = {1.0,2.0,1.0,4.0}; + boost::random::dynamic_discrete_distribution<> dist(probs); + boost::random::dynamic_discrete_distribution<>::param_type param = dist.param(); + CHECK_PROBABILITIES(param.probabilities(), (std::vector{.125,.25,.125,.5})); + boost::random::dynamic_discrete_distribution<> copy1(param); + BOOST_CHECK_EQUAL(dist, copy1); + boost::random::dynamic_discrete_distribution<> copy2; + copy2.param(param); + BOOST_CHECK_EQUAL(dist, copy2); + + boost::random::dynamic_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::dynamic_discrete_distribution<>::param_type param_default; + CHECK_PROBABILITIES(param_default.probabilities(), (std::vector{1.0})); + BOOST_CHECK(param != param_default); + BOOST_CHECK(!(param == param_default)); + +#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST + boost::random::dynamic_discrete_distribution<>::param_type + parm_il = (std::vector{ 1, 2, 1, 4 }); + CHECK_PROBABILITIES(parm_il.probabilities(), (std::vector{.125,.25,.125,.5})); +#endif + + boost::random::dynamic_discrete_distribution<>::param_type parm_r(probs); + CHECK_PROBABILITIES(parm_r.probabilities(), (std::vector{.125,.25,.125,.5})); + + boost::random::dynamic_discrete_distribution<>::param_type + parm_it(probs.begin(), probs.end()); + CHECK_PROBABILITIES(parm_it.probabilities(), (std::vector{.125,.25,.125,.5})); + + boost::random::dynamic_discrete_distribution<>::param_type + parm_fun(4, 99, 115, gen()); + CHECK_PROBABILITIES(parm_fun.probabilities(), (std::vector{.125,.25,.125,.5})); + +} + +BOOST_AUTO_TEST_CASE(test_min_max) { + std::vector probs = {1.0,2.0,1.0}; + boost::random::dynamic_discrete_distribution<> dist; + BOOST_CHECK_EQUAL((dist.min)(), 0); + BOOST_CHECK_EQUAL((dist.max)(), 0); + boost::random::dynamic_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 probs = {1.0,2.0,1.0,4.0}; + boost::random::dynamic_discrete_distribution<> dist; + boost::random::dynamic_discrete_distribution<> dist_copy(dist); + boost::random::dynamic_discrete_distribution<> dist_r(probs); + boost::random::dynamic_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 probs = {1.0,2.0,1.0,4.0}; + boost::random::dynamic_discrete_distribution<> dist(probs); + std::stringstream stream; + stream << dist; + boost::random::dynamic_discrete_distribution<> restored_dist; + stream >> restored_dist; + BOOST_CHECK_EQUAL(dist, restored_dist); +} + +BOOST_AUTO_TEST_CASE(test_generation) { + std::vector probs = {0.0,1.0}; + std::minstd_rand0 gen; + boost::random::dynamic_discrete_distribution<> dist; + boost::random::dynamic_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); + } +} + +template +double testSelection( int n, int max){ + std::uniform_real_distribution uniform_real_dist(.5,1.0); + + max = std::min(max,10*n); + boost::mt19937 generator; + std::vector weights(n); + double sum = 0; + for (int i=0;i dist(weights); + std::vector results(n,0); + for (int i=0;i +double testUpdate( int n, int max){ + std::uniform_real_distribution uniform_real_dist(.5,1.0); + + max = std::min(max,10*n); + boost::mt19937 generator; + std::vector weights(n); + double sum = 0; + for (int i=0;i dist(weights); + std::uniform_int_distribution uniform_dist(0,n-1); + for (int i=0;i results(n,0); + for (int i=0;i +// bool stat_test(boost::random::dynamic_discrete_distribution distrubiton, std::vector weights){ +// boost::mt19937 generator; +// auto probabilities = weights; +// int sum = 0; +// for (int i=0; i dist({0,1,2,3,4,5}); + BOOST_CHECK(dist.get_weight(0)==0); + BOOST_CHECK(dist.get_weight(5)==5); + +} + +//Question should we throw errors when container is empty? +//Question, when container is empty, should we throw errors during selection? +BOOST_AUTO_TEST_CASE(adding_removing){ + boost::random::dynamic_discrete_distribution dist{0,1,2,3,4}; + dist.pop_back(5); + BOOST_CHECK(dist.total_weight()==0); + + std::uniform_real_distribution uniform_real_dist(.5,1.0); + boost::mt19937 generator; + std::vector weights(35); + double sum = 0; + for (int i=0;i<35;i++){ + double rand = uniform_real_dist(generator); + dist.push_back(rand); + weights[i] = rand; + sum+=rand; + } + for (size_t i=0;i results(35); + int iters = 1000; + for (int i=0;i newResults(weights.size()); + for (int i=0; i(2,1000)<0.99); + BOOST_CHECK(testSelection<2>(3,1000)<0.99); + BOOST_CHECK(testSelection<2>(30,100000)<0.99); + //k=8 + BOOST_CHECK(testSelection<8>(2,1000)<0.99); + BOOST_CHECK(testSelection<8>(9,1000)<0.99); + BOOST_CHECK(testSelection<8>(100,100000)<0.99); + //k= 32 + BOOST_CHECK(testSelection<32>(2,1000)<0.99); + BOOST_CHECK(testSelection<32>(33,1000)<0.99); + BOOST_CHECK(testSelection<32>(1000,1000000)<0.99); + +} + +BOOST_AUTO_TEST_CASE(update){ + //k=2 + BOOST_CHECK(testUpdate<2>(2,1000)<0.99); + BOOST_CHECK(testUpdate<2>(3,1000)<0.99); + BOOST_CHECK(testUpdate<2>(30,100000)<0.99); + //k=8 + BOOST_CHECK(testUpdate<8>(2,1000)<0.99); + BOOST_CHECK(testUpdate<8>(9,1000)<0.99); + BOOST_CHECK(testUpdate<8>(100,10000)<0.99); + //k= 32 + BOOST_CHECK(testUpdate<32>(2,1000)<0.99); + BOOST_CHECK(testUpdate<32>(33,100000)<0.99); + BOOST_CHECK(testUpdate<32>(1000,100000)<0.99); + +}