Chaste  Release::2018.1
Boost165GammaDistribution.hpp
Go to the documentation of this file.
1 
17 /* boost random/gamma_distribution.hpp header file
18  *
19  * Copyright Jens Maurer 2002
20  * Copyright Steven Watanabe 2010
21  * Distributed under the Boost Software License, Version 1.0. (See
22  * accompanying file LICENSE_1_0.txt or copy at
23  * http://www.boost.org/LICENSE_1_0.txt)
24  *
25  * See http://www.boost.org for most recent version including documentation.
26  *
27  * $Id$
28  *
29  */
30 
31 #ifndef BOOST_165_RANDOM_GAMMA_DISTRIBUTION_HPP
32 #define BOOST_165_RANDOM_GAMMA_DISTRIBUTION_HPP
33 
34 #include <boost/assert.hpp>
35 #include <boost/config/no_tr1/cmath.hpp>
36 #include <boost/limits.hpp>
37 #include <boost/random/detail/config.hpp>
38 //#include <boost/random/exponential_distribution.hpp> // swapped for chaste copy of exponential v1.65
39 #include <boost/static_assert.hpp>
40 #include <iosfwd>
41 #include <istream>
43 
44 namespace boost
45 {
46 namespace random
47 {
48 
49  // The algorithm is taken from Knuth
50 
58  template <class RealType = double>
60  {
61  public:
62  typedef RealType input_type;
63  typedef RealType result_type;
64 
65  class param_type
66  {
67  public:
69 
76  param_type(const RealType& alpha_arg = RealType(1.0),
77  const RealType& beta_arg = RealType(1.0))
78  : _alpha(alpha_arg), _beta(beta_arg)
79  {
80  }
81 
83  RealType alpha() const { return _alpha; }
85  RealType beta() const { return _beta; }
86 
87 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
88 
89  template <class CharT, class Traits>
90  friend std::basic_ostream<CharT, Traits>&
91  operator<<(std::basic_ostream<CharT, Traits>& os,
92  const param_type& parm)
93  {
94  os << parm._alpha << ' ' << parm._beta;
95  return os;
96  }
97 
99  template <class CharT, class Traits>
100  friend std::basic_istream<CharT, Traits>&
101  operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm)
102  {
103  is >> parm._alpha >> std::ws >> parm._beta;
104  return is;
105  }
106 #endif
107 
109  friend bool operator==(const param_type& lhs, const param_type& rhs)
110  {
111  return lhs._alpha == rhs._alpha && lhs._beta == rhs._beta;
112  }
114  friend bool operator!=(const param_type& lhs, const param_type& rhs)
115  {
116  return !(lhs == rhs);
117  }
118 
119  private:
120  RealType _alpha;
121  RealType _beta;
122  };
123 
124 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
125  BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
126 #endif
127 
133  explicit gamma_distribution_v165(const result_type& alpha_arg = result_type(1.0),
134  const result_type& beta_arg = result_type(1.0))
135  : _exp(result_type(1)), _alpha(alpha_arg), _beta(beta_arg)
136  {
137  BOOST_ASSERT(_alpha > result_type(0));
138  BOOST_ASSERT(_beta > result_type(0));
139  init();
140  }
141 
143  explicit gamma_distribution_v165(const param_type& parm)
144  : _exp(result_type(1)), _alpha(parm.alpha()), _beta(parm.beta())
145  {
146  init();
147  }
148 
149  // compiler-generated copy ctor and assignment operator are fine
150 
152  RealType alpha() const { return _alpha; }
154  RealType beta() const { return _beta; }
156  RealType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 0; }
157  /* Returns the largest value that the distribution can produce. */
158  RealType max BOOST_PREVENT_MACRO_SUBSTITUTION() const
159  {
160  return (std::numeric_limits<RealType>::infinity)();
161  }
162 
164  param_type param() const { return param_type(_alpha, _beta); }
166  void param(const param_type& parm)
167  {
168  _alpha = parm.alpha();
169  _beta = parm.beta();
170  init();
171  }
172 
177  void reset() { _exp.reset(); }
178 
183  template <class Engine>
184  result_type operator()(Engine& eng)
185  {
186 #ifndef BOOST_NO_STDC_NAMESPACE
187  // allow for Koenig lookup
188  using std::tan;
189  using std::sqrt;
190  using std::exp;
191  using std::log;
192  using std::pow;
193 #endif
194  if (_alpha == result_type(1))
195  {
196  return _exp(eng) * _beta;
197  }
198  else if (_alpha > result_type(1))
199  {
200  // Can we have a boost::mathconst please?
201  const result_type pi = result_type(3.14159265358979323846);
202  for (;;)
203  {
204  result_type y = tan(pi * uniform_01<RealType>()(eng));
205  result_type x = sqrt(result_type(2) * _alpha - result_type(1)) * y
206  + _alpha - result_type(1);
207  if (x <= result_type(0))
208  continue;
209  if (uniform_01<RealType>()(eng) > (result_type(1) + y * y) * exp((_alpha - result_type(1))
210  * log(x / (_alpha - result_type(1)))
211  - sqrt(result_type(2) * _alpha
212  - result_type(1))
213  * y))
214  continue;
215  return x * _beta;
216  }
217  }
218  else /* alpha < 1.0 */
219  {
220  for (;;)
221  {
222  result_type u = uniform_01<RealType>()(eng);
223  result_type y = _exp(eng);
224  result_type x, q;
225  if (u < _p)
226  {
227  x = exp(-y / _alpha);
228  q = _p * exp(-x);
229  }
230  else
231  {
232  x = result_type(1) + y;
233  q = _p + (result_type(1) - _p) * pow(x, _alpha - result_type(1));
234  }
235  if (u >= q)
236  continue;
237  return x * _beta;
238  }
239  }
240  }
241 
242  template <class URNG>
243  RealType operator()(URNG& urng, const param_type& parm) const
244  {
245  return gamma_distribution_v165(parm)(urng);
246  }
247 
248 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
249 
250  template <class CharT, class Traits>
251  friend std::basic_ostream<CharT, Traits>&
252  operator<<(std::basic_ostream<CharT, Traits>& os,
253  const gamma_distribution_v165& gd)
254  {
255  os << gd.param();
256  return os;
257  }
258 
260  template <class CharT, class Traits>
261  friend std::basic_istream<CharT, Traits>&
262  operator>>(std::basic_istream<CharT, Traits>& is, gamma_distribution_v165& gd)
263  {
264  gd.read(is);
265  return is;
266  }
267 #endif
268 
273  friend bool operator==(const gamma_distribution_v165& lhs,
274  const gamma_distribution_v165& rhs)
275  {
276  return lhs._alpha == rhs._alpha
277  && lhs._beta == rhs._beta
278  && lhs._exp == rhs._exp;
279  }
280 
285  friend bool operator!=(const gamma_distribution_v165& lhs,
286  const gamma_distribution_v165& rhs)
287  {
288  return !(lhs == rhs);
289  }
290 
291  private:
293 
294  template <class CharT, class Traits>
295  void read(std::basic_istream<CharT, Traits>& is)
296  {
297  param_type parm;
298  if (is >> parm)
299  {
300  param(parm);
301  }
302  }
303 
304  void init()
305  {
306 #ifndef BOOST_NO_STDC_NAMESPACE
307  // allow for Koenig lookup
308  using std::exp;
309 #endif
310  _p = exp(result_type(1)) / (_alpha + exp(result_type(1)));
311  }
313 
314  exponential_distribution_v165<RealType> _exp;
315  result_type _alpha;
316  result_type _beta;
317  // some data precomputed from the parameters
318  result_type _p;
319  };
320 
321 } // namespace random
322 
323 } // namespace boost
324 
325 #endif // BOOST_165_RANDOM_GAMMA_DISTRIBUTION_HPP
friend std::basic_istream< CharT, Traits > & operator>>(std::basic_istream< CharT, Traits > &is, param_type &parm)
friend bool operator==(const gamma_distribution_v165 &lhs, const gamma_distribution_v165 &rhs)
param_type(const RealType &alpha_arg=RealType(1.0), const RealType &beta_arg=RealType(1.0))
gamma_distribution_v165(const result_type &alpha_arg=result_type(1.0), const result_type &beta_arg=result_type(1.0))
friend bool operator!=(const gamma_distribution_v165 &lhs, const gamma_distribution_v165 &rhs)
friend bool operator!=(const param_type &lhs, const param_type &rhs)
friend bool operator==(const param_type &lhs, const param_type &rhs)
friend std::basic_istream< CharT, Traits > & operator>>(std::basic_istream< CharT, Traits > &is, gamma_distribution_v165 &gd)