9 #ifndef IMPKERNEL_PIECEWISE_LINEAR_DISTRIBUTION_H
10 #define IMPKERNEL_PIECEWISE_LINEAR_DISTRIBUTION_H
12 #include <IMP/kernel_config.h>
14 #include <boost/version.hpp>
15 #include <boost/random/uniform_real.hpp>
23 IMPKERNEL_BEGIN_NAMESPACE
34 template <
class T =
double>
39 mutable boost::uniform_real<> un_;
40 double get_divider(
int divider)
const {
return dividers_[divider]; }
41 double get_accum(
int divider)
const {
return accum_[divider]; }
42 double get_slope(
int divider)
const {
43 return (weights_[divider + 1] - weights_[divider]) /
44 (dividers_[divider + 1] - dividers_[divider]);
46 double get_weight(
int divider)
const {
return weights_[divider]; }
51 dividers_.push_back(0);
52 dividers_.push_back(1);
55 weights_.push_back(1);
56 weights_.push_back(1);
63 template <
class LIt,
class WIt>
66 : dividers_(locations_begin, locations_end),
67 accum_(std::distance(locations_begin, locations_end)),
68 weights_(std::distance(locations_begin, locations_end)),
70 for (
unsigned int i = 0; i < weights_.size(); ++i) {
71 weights_[i] = *weights_begin;
75 for (
unsigned int i = 1; i < accum_.size(); ++i) {
76 double wid = (dividers_[i] - dividers_[i - 1]);
77 accum_[i] = accum_[i - 1] + wid * .5 * (weights_[i - 1] + weights_[i]);
79 double total = accum_.back();
80 #if IMP_HAS_CHECKS >= IMP_INTERNAL
81 for (
unsigned int i = 1; i < dividers_.size(); ++i) {
82 if (dividers_[i] <= dividers_[i - 1]) {
83 std::cerr <<
"Found non-monotonic locations: ";
84 std::copy(locations_begin, locations_end,
85 std::ostream_iterator<double>(std::cerr,
", "));
86 std::cerr << std::endl;
87 throw std::logic_error(
"Locations are not strictly increasing");
91 if (std::abs(total - 1.0) > .05) {
92 std::cerr <<
"Found weight " << total << std::endl;
93 std::cerr <<
"Locations: [";
94 std::copy(locations_begin, locations_end,
95 std::ostream_iterator<double>(std::cerr,
", "));
96 std::cerr <<
"]" << std::endl;
97 std::cerr <<
"Weights: [";
98 std::copy(weights_.begin(), weights_.end(),
99 std::ostream_iterator<double>(std::cerr,
", "));
100 std::cerr <<
"]" << std::endl;
101 throw std::logic_error(
102 "The total weight of the distribution is not close to 1");
106 for (
unsigned int i = 0; i < accum_.size(); ++i) {
122 double operator()(RNG& rng)
const {
123 double rn = un_(rng);
125 std::upper_bound(accum_.begin() + 1, accum_.end(), rn);
126 if (it == accum_.end()) {
127 std::cerr <<
"Hit end " << rn <<
" " << accum_.back() << std::endl;
128 return dividers_.back();
130 int offset = it - accum_.begin() - 1;
131 double pd = get_weight(offset);
132 double ad = get_accum(offset);
133 double sd = get_slope(offset);
134 double leftd = rn - ad;
135 double radix = std::sqrt(pd * pd + 2 * leftd * sd);
138 xmd = (-pd + radix) / sd;
140 xmd = (-pd + radix) / (sd);
144 std::cerr <<
"negative " << std::endl;
149 return xmd + get_divider(offset);
154 IMPKERNEL_END_NAMESPACE
piecewise_linear_distribution(LIt locations_begin, LIt locations_end, WIt weights_begin)
piecewise_linear_distribution()
construct a uniform 01
A more IMP-like version of the std::vector.
A class for storing lists of IMP items.