IMP  2.0.0
The Integrative Modeling Platform
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-2013 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 
25 /** Draw random numbers from a distribution defined as a piecewise
26  linear function. It models the boost random number generators
27  and is made redundant by the boost::piecewise_linear_distribution in
28  boost 1.47. Currently, it won't use that class if available, but
29  that is just due to lack of boost 1.47 to test on.
30 */
31 template <class T=double>
33  Vector<T> dividers_;
34  Vector<T> accum_;
35  Vector<T> weights_;
36  mutable boost::uniform_real<> un_;
37  double get_divider(int divider) const {
38  return dividers_[divider];
39  }
40  double get_accum(int divider) const {
41  return accum_[divider];
42  }
43  double get_slope(int divider) const {
44  return (weights_[divider+1]-weights_[divider])
45  /(dividers_[divider+1]-dividers_[divider]);
46  }
47  double get_weight(int divider) const {
48  return weights_[divider]; }
49  public:
50  //! construct a uniform 01
52  dividers_.push_back(0);
53  dividers_.push_back(1);
54  accum_.push_back(0);
55  accum_.push_back(1);
56  weights_.push_back(1);
57  weights_.push_back(1);
58  }
59  /** Construct the distribution by interpolating between samples
60  at locations[i] with corresponding weight weights[i]. The
61  weight outside of the first and last location is assume to
62  be immediately 0.
63  */
64  template <class LIt, class WIt>
65  piecewise_linear_distribution(LIt locations_begin, LIt locations_end,
66  WIt weights_begin) :
67  dividers_(locations_begin,
68  locations_end),
69  accum_(std::distance(locations_begin, locations_end)),
70  weights_(std::distance(locations_begin, locations_end)),
71  un_(0,1) {
72  for (unsigned int i=0; i< weights_.size(); ++i) {
73  weights_[i]=*weights_begin;
74  ++weights_begin;
75  }
76  accum_[0]=0;
77  for (unsigned int i=1; i< accum_.size(); ++i) {
78  double wid=(dividers_[i]-dividers_[i-1]);
79  accum_[i]= accum_[i-1]+wid*.5*(weights_[i-1]+weights_[i]);
80  }
81  double total=accum_.back();
82 #if IMP_HAS_CHECKS >= IMP_INTERNAL
83  for (unsigned int i=1; i< dividers_.size(); ++i) {
84  if (dividers_[i] <= dividers_[i-1]) {
85  std::cerr << "Found non-monotonic locations: ";
86  std::copy(locations_begin, locations_end,
87  std::ostream_iterator<double>(std::cerr, ", "));
88  std::cerr << std::endl;
89  throw std::logic_error("Locations are not strictly increasing");
90  }
91  }
92  // check that the total weight is ~1
93  if (std::abs(total-1.0) > .05){
94  std::cerr << "Found weight " << total << std::endl;
95  std::cerr << "Locations: [";
96  std::copy(locations_begin, locations_end,
97  std::ostream_iterator<double>(std::cerr, ", "));
98  std::cerr << "]" << std::endl;
99  std::cerr << "Weights: [";
100  std::copy(weights_.begin(), weights_.end(),
101  std::ostream_iterator<double>(std::cerr, ", "));
102  std::cerr << "]" << std::endl;
103  throw
104  std::logic_error
105  ("The total weight of the distribution is not close to 1");
106  }
107  #endif
108  // correct for discretization errors
109  for (unsigned int i=0; i< accum_.size(); ++i) {
110  accum_[i]/=total;
111  }
112  /*std::cout << "dividers: ";
113  std::copy(dividers_.begin(), dividers_.end(),
114  std::ostream_iterator<double>(std::cout, ", "));
115  std::cout << std::endl << "sums: ";
116  std::copy(accum_.begin(), accum_.end(),
117  std::ostream_iterator<double>(std::cout, ", "));
118  std::cout << std::endl;
119  std::cout << "weights: ";
120  std::copy(weights_.begin(), weights_.end(),
121  std::ostream_iterator<double>(std::cout, ", "));
122  std::cout << std::endl;*/
123  }
124  template <class RNG>
125  double operator()(RNG& rng) const {
126  double rn= un_(rng);
127  typename Vector<T>::const_iterator it
128  = std::upper_bound(accum_.begin()+1, accum_.end(), rn);
129  if (it == accum_.end()) {
130  std::cerr << "Hit end " << rn << " " << accum_.back() << std::endl;
131  return dividers_.back();
132  } else {
133  int offset= it-accum_.begin()-1;
134  double pd= get_weight(offset);
135  double ad= get_accum(offset);
136  double sd= get_slope(offset);
137  double leftd=rn-ad;
138  double radix=std::sqrt(pd*pd+2*leftd*sd);
139  double xmd;
140  if (sd < 0) {
141  xmd= (-pd + radix)/sd;
142  } else {
143  xmd= (-pd + radix)/(sd);
144  }
145 #ifndef NDEBUG
146  if (xmd <0) {
147  std::cerr << "negative " << std::endl;
148  }
149 #endif
150  /*std::cout << rn << " " << pd << " " << ad << " " << sd << " "
151  << leftd << " " << xmd << std::endl;*/
152  return xmd+get_divider(offset);
153  }
154  }
155 };
156 
157 
158 
159 IMPBASE_END_NAMESPACE
160 
161 #endif /* IMPBASE_PIECEWISE_LINEAR_DISTRIBUTION_H */