IMP  2.2.1
The Integrative Modeling Platform
base/piecewise_linear_distribution.h
Go to the documentation of this file.
1 /**
2  * \file IMP/base/piecewise_linear_distribution.h
3  * \brief boost piecewise linear.
4  *
5  * Copyright 2007-2014 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPBASE_PIECEWISE_LINEAR_DISTRIBUTION_H
10 #define IMPBASE_PIECEWISE_LINEAR_DISTRIBUTION_H
11 
12 #include <IMP/base/base_config.h>
13 #include "Vector.h"
14 #include <boost/version.hpp>
15 #include <boost/random/uniform_real.hpp>
16 #include <numeric>
17 #include <algorithm>
18 #include <cmath>
19 #include <iostream>
20 #include <iterator>
21 #include <stdexcept>
22 
23 IMPBASE_BEGIN_NAMESPACE
24 /** Draw random numbers from a distribution defined as a
25  piecewise
26  linear function. It models the boost random number
27  generators
28  and is made redundant by the
29  boost::piecewise_linear_distribution in
30  boost 1.47. Currently, it won't use that class if
31  available, but
32  that is just due to lack of boost 1.47 to test on.
33  */
34 template <class T = double>
36  Vector<T> dividers_;
37  Vector<T> accum_;
38  Vector<T> weights_;
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]);
45  }
46  double get_weight(int divider) const { return weights_[divider]; }
47 
48  public:
49  //! construct a uniform 01
51  dividers_.push_back(0);
52  dividers_.push_back(1);
53  accum_.push_back(0);
54  accum_.push_back(1);
55  weights_.push_back(1);
56  weights_.push_back(1);
57  }
58  /** Construct the distribution by interpolating between samples
59  at locations[i] with corresponding weight weights[i]. The
60  weight outside of the first and last location is assume to
61  be immediately 0.
62  */
63  template <class LIt, class WIt>
64  piecewise_linear_distribution(LIt locations_begin, LIt locations_end,
65  WIt weights_begin)
66  : dividers_(locations_begin, locations_end),
67  accum_(std::distance(locations_begin, locations_end)),
68  weights_(std::distance(locations_begin, locations_end)),
69  un_(0, 1) {
70  for (unsigned int i = 0; i < weights_.size(); ++i) {
71  weights_[i] = *weights_begin;
72  ++weights_begin;
73  }
74  accum_[0] = 0;
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]);
78  }
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");
88  }
89  }
90  // check that the total weight is ~1
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");
103  }
104 #endif
105  // correct for discretization errors
106  for (unsigned int i = 0; i < accum_.size(); ++i) {
107  accum_[i] /= total;
108  }
109  /*std::cout << "dividers: ";
110  std::copy(dividers_.begin(), dividers_.end(),
111  std::ostream_iterator<double>(std::cout, ", "));
112  std::cout << std::endl << "sums: ";
113  std::copy(accum_.begin(), accum_.end(),
114  std::ostream_iterator<double>(std::cout, ", "));
115  std::cout << std::endl;
116  std::cout << "weights: ";
117  std::copy(weights_.begin(), weights_.end(),
118  std::ostream_iterator<double>(std::cout, ", "));
119  std::cout << std::endl;*/
120  }
121  template <class RNG>
122  double operator()(RNG& rng) const {
123  double rn = un_(rng);
124  typename Vector<T>::const_iterator it =
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();
129  } else {
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);
136  double xmd;
137  if (sd < 0) {
138  xmd = (-pd + radix) / sd;
139  } else {
140  xmd = (-pd + radix) / (sd);
141  }
142 #ifndef NDEBUG
143  if (xmd < 0) {
144  std::cerr << "negative " << std::endl;
145  }
146 #endif
147  /*std::cout << rn << " " << pd << " " << ad << " " << sd << " "
148  << leftd << " " << xmd << std::endl;*/
149  return xmd + get_divider(offset);
150  }
151  }
152 };
153 
154 IMPBASE_END_NAMESPACE
155 
156 #endif /* IMPBASE_PIECEWISE_LINEAR_DISTRIBUTION_H */
piecewise_linear_distribution(LIt locations_begin, LIt locations_end, WIt weights_begin)
A class for storing lists of IMP items.