Chaste  Release::2018.1
Boost165ExponentialDistribution.hpp
Go to the documentation of this file.
1 
17 /* boost random/exponential_distribution.hpp header file
18  *
19  * Copyright Jens Maurer 2000-2001
20  * Copyright Steven Watanabe 2011
21  * Copyright Jason Rhinelander 2016
22  * Distributed under the Boost Software License, Version 1.0. (See
23  * accompanying file LICENSE_1_0.txt or copy at
24  * http://www.boost.org/LICENSE_1_0.txt)
25  *
26  * See http://www.boost.org for most recent version including documentation.
27  *
28  * $Id$
29  *
30  * Revision history
31  * 2001-02-18 moved to individual header files
32  */
33 
34 #ifndef BOOST_165_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP
35 #define BOOST_165_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP
36 
37 #include <boost/assert.hpp>
38 #include <boost/config/no_tr1/cmath.hpp>
39 #include <boost/limits.hpp>
40 #include <boost/random/detail/config.hpp>
41 // #include <boost/random/detail/int_float_pair.hpp // Not in older boosts - copied into Chaste.
42 #include <boost/random/detail/operators.hpp>
43 #include <boost/random/uniform_01.hpp>
44 #include <iosfwd>
45 
46 #include "Boost165IntFloatPair.hpp" // version 1.65 of int_float_pair.hpp packaged with Chaste.
47 
48 namespace boost
49 {
50 namespace random
51 {
52 
53  namespace detail
54  {
55 
56  // tables for the ziggurat algorithm
57  template <class RealType>
59  {
60  static const RealType table_x[257];
61  static const RealType table_y[257];
62  };
63 
64  template <class RealType>
65  const RealType exponential_table<RealType>::table_x[257] = {
66  8.6971174701310497140, 7.6971174701310497140, 6.9410336293772123602, 6.4783784938325698538,
67  6.1441646657724730491, 5.8821443157953997963, 5.6664101674540337371, 5.4828906275260628694,
68  5.3230905057543986131, 5.1814872813015010392, 5.0542884899813047117, 4.9387770859012514838,
69  4.8329397410251125881, 4.7352429966017412526, 4.6444918854200854873, 4.5597370617073515513,
70  4.4802117465284221949, 4.4052876934735729805, 4.3344436803172730116, 4.2672424802773661873,
71  4.2033137137351843802, 4.1423408656640511251, 4.0840513104082974638, 4.0282085446479365106,
72  3.9746060666737884793, 3.9230625001354895926, 3.8734176703995089983, 3.8255294185223367372,
73  3.7792709924116678992, 3.7345288940397975350, 3.6912010902374189454, 3.6491955157608538478,
74  3.6084288131289096339, 3.5688252656483374051, 3.5303158891293438633, 3.4928376547740601814,
75  3.4563328211327607625, 3.4207483572511205323, 3.3860354424603017887, 3.3521490309001100106,
76  3.3190474709707487166, 3.2866921715990692095, 3.2550473085704501813, 3.2240795652862645207,
77  3.1937579032122407483, 3.1640533580259734580, 3.1349388580844407393, 3.1063890623398246660,
78  3.0783802152540905188, 3.0508900166154554479, 3.0238975044556767713, 2.9973829495161306949,
79  2.9713277599210896472, 2.9457143948950456386, 2.9205262865127406647, 2.8957477686001416838,
80  2.8713640120155362592, 2.8473609656351888266, 2.8237253024500354905, 2.8004443702507381944,
81  2.7775061464397572041, 2.7548991965623453650, 2.7326126361947007411, 2.7106360958679293686,
82  2.6889596887418041593, 2.6675739807732670816, 2.6464699631518093905, 2.6256390267977886123,
83  2.6050729387408355373, 2.5847638202141406911, 2.5647041263169053687, 2.5448866271118700928,
84  2.5253043900378279427, 2.5059507635285939648, 2.4868193617402096807, 2.4679040502973649846,
85  2.4491989329782498908, 2.4306983392644199088, 2.4123968126888708336, 2.3942890999214583288,
86  2.3763701405361408194, 2.3586350574093374601, 2.3410791477030346875, 2.3236978743901964559,
87  2.3064868582835798692, 2.2894418705322694265, 2.2725588255531546952, 2.2558337743672190441,
88  2.2392628983129087111, 2.2228425031110364013, 2.2065690132576635755, 2.1904389667232199235,
89  2.1744490099377744673, 2.1585958930438856781, 2.1428764653998416425, 2.1272876713173679737,
90  2.1118265460190418108, 2.0964902118017147637, 2.0812758743932248696, 2.0661808194905755036,
91  2.0512024094685848641, 2.0363380802487695916, 2.0215853383189260770, 2.0069417578945183144,
92  1.9924049782135764992, 1.9779727009573602295, 1.9636426877895480401, 1.9494127580071845659,
93  1.9352807862970511135, 1.9212447005915276767, 1.9073024800183871196, 1.8934521529393077332,
94  1.8796917950722108462, 1.8660195276928275962, 1.8524335159111751661, 1.8389319670188793980,
95  1.8255131289035192212, 1.8121752885263901413, 1.7989167704602903934, 1.7857359354841254047,
96  1.7726311792313049959, 1.7596009308890742369, 1.7466436519460739352, 1.7337578349855711926,
97  1.7209420025219350428, 1.7081947058780575683, 1.6955145241015377061, 1.6829000629175537544,
98  1.6703499537164519163, 1.6578628525741725325, 1.6454374393037234057, 1.6330724165359912048,
99  1.6207665088282577216, 1.6085184617988580769, 1.5963270412864831349, 1.5841910325326886695,
100  1.5721092393862294810, 1.5600804835278879161, 1.5481036037145133070, 1.5361774550410318943,
101  1.5243009082192260050, 1.5124728488721167573, 1.5006921768428164936, 1.4889578055167456003,
102  1.4772686611561334579, 1.4656236822457450411, 1.4540218188487932264, 1.4424620319720121876,
103  1.4309432929388794104, 1.4194645827699828254, 1.4080248915695353509, 1.3966232179170417110,
104  1.3852585682631217189, 1.3739299563284902176, 1.3626364025050864742, 1.3513769332583349176,
105  1.3401505805295045843, 1.3289563811371163220, 1.3177933761763245480, 1.3066606104151739482,
106  1.2955571316866007210, 1.2844819902750125450, 1.2734342382962410994, 1.2624129290696153434,
107  1.2514171164808525098, 1.2404458543344064544, 1.2294981956938491599, 1.2185731922087903071,
108  1.2076698934267612830, 1.1967873460884031665, 1.1859245934042023557, 1.1750806743109117687,
109  1.1642546227056790397, 1.1534454666557748056, 1.1426522275816728928, 1.1318739194110786733,
110  1.1211095477013306083, 1.1103581087274114281, 1.0996185885325976575, 1.0888899619385472598,
111  1.0781711915113727024, 1.0674612264799681530, 1.0567590016025518414, 1.0460634359770445503,
112  1.0353734317905289496, 1.0246878730026178052, 1.0140056239570971074, 1.0033255279156973717,
113  0.99264640550727647009, 0.98196705308506317914, 0.97128624098390397896, 0.96060271166866709917,
114  0.94991517776407659940, 0.93922231995526297952, 0.92852278474721113999, 0.91781518207004493915,
115  0.90709808271569100600, 0.89637001558989069006, 0.88562946476175228052, 0.87487486629102585352,
116  0.86410460481100519511, 0.85331700984237406386, 0.84251035181036928333, 0.83168283773427388393,
117  0.82083260655441252290, 0.80995772405741906620, 0.79905617735548788109, 0.78812586886949324977,
118  0.77716460975913043936, 0.76617011273543541328, 0.75513998418198289808, 0.74407171550050873971,
119  0.73296267358436604916, 0.72181009030875689912, 0.71061105090965570413, 0.69936248110323266174,
120  0.68806113277374858613, 0.67670356802952337911, 0.66528614139267855405, 0.65380497984766565353,
121  0.64225596042453703448, 0.63063468493349100113, 0.61893645139487678178, 0.60715622162030085137,
122  0.59528858429150359384, 0.58332771274877027785, 0.57126731653258903915, 0.55910058551154127652,
123  0.54682012516331112550, 0.53441788123716615385, 0.52188505159213564105, 0.50921198244365495319,
124  0.49638804551867159754, 0.48340149165346224782, 0.47023927508216945338, 0.45688684093142071279,
125  0.44332786607355296305, 0.42954394022541129589, 0.41551416960035700100, 0.40121467889627836229,
126  0.38661797794112021568, 0.37169214532991786118, 0.35639976025839443721, 0.34069648106484979674,
127  0.32452911701691008547, 0.30783295467493287307, 0.29052795549123115167, 0.27251318547846547924,
128  0.25365836338591284433, 0.23379048305967553619, 0.21267151063096745264, 0.18995868962243277774,
129  0.16512762256418831796, 0.13730498094001380420, 0.10483850756582017915, 0.063852163815003480173,
130  0
131  };
132 
133  template <class RealType>
134  const RealType exponential_table<RealType>::table_y[257] = {
135  0, 0.00045413435384149675545, 0.00096726928232717452884, 0.0015362997803015723824,
136  0.0021459677437189061793, 0.0027887987935740759640, 0.0034602647778369039855, 0.0041572951208337952532,
137  0.0048776559835423925804, 0.0056196422072054831710, 0.0063819059373191794422, 0.0071633531836349841425,
138  0.0079630774380170392396, 0.0087803149858089752347, 0.0096144136425022094101, 0.010464810181029979488,
139  0.011331013597834597488, 0.012212592426255380661, 0.013109164931254991070, 0.014020391403181937334,
140  0.014945968011691148079, 0.015885621839973162490, 0.016839106826039946359, 0.017806200410911360563,
141  0.018786700744696029497, 0.019780424338009741737, 0.020787204072578117603, 0.021806887504283582125,
142  0.022839335406385238829, 0.023884420511558170348, 0.024942026419731782971, 0.026012046645134218076,
143  0.027094383780955798424, 0.028188948763978634421, 0.029295660224637394015, 0.030414443910466605492,
144  0.031545232172893605499, 0.032687963508959533317, 0.033842582150874329031, 0.035009037697397411067,
145  0.036187284781931419754, 0.037377282772959360128, 0.038578995503074859626, 0.039792391023374122670,
146  0.041017441380414820816, 0.042254122413316231413, 0.043502413568888183301, 0.044762297732943280694,
147  0.046033761076175166762, 0.047316792913181548703, 0.048611385573379494401, 0.049917534282706374944,
148  0.051235237055126279830, 0.052564494593071689595, 0.053905310196046085104, 0.055257689676697038322,
149  0.056621641283742874438, 0.057997175631200659098, 0.059384305633420264487, 0.060783046445479636051,
150  0.062193415408540996150, 0.063615431999807331076, 0.065049117786753755036, 0.066494496385339779043,
151  0.067951593421936607770, 0.069420436498728751675, 0.070901055162371828426, 0.072393480875708743023,
152  0.073897746992364746308, 0.075413888734058408453, 0.076941943170480510100, 0.078481949201606426042,
153  0.080033947542319910023, 0.081597980709237420930, 0.083174093009632380354, 0.084762330532368125386,
154  0.086362741140756912277, 0.087975374467270219300, 0.089600281910032864534, 0.091237516631040162057,
155  0.092887133556043546523, 0.094549189376055853718, 0.096223742550432800103, 0.097910853311492199618,
156  0.099610583670637128826, 0.10132299742595363588, 0.10304816017125771553, 0.10478613930657016928,
157  0.10653700405000166218, 0.10830082545103379867, 0.11007767640518539026, 0.11186763167005629731,
158  0.11367076788274431301, 0.11548716357863353664, 0.11731689921155557057, 0.11916005717532768467,
159  0.12101672182667483729, 0.12288697950954513498, 0.12477091858083096578, 0.12666862943751066518,
160  0.12858020454522817870, 0.13050573846833078225, 0.13244532790138752023, 0.13439907170221363078,
161  0.13636707092642885841, 0.13834942886358021406, 0.14034625107486244210, 0.14235764543247220043,
162  0.14438372216063476473, 0.14642459387834493787, 0.14848037564386679222, 0.15055118500103990354,
163  0.15263714202744286154, 0.15473836938446807312, 0.15685499236936522013, 0.15898713896931420572,
164  0.16113493991759203183, 0.16329852875190180795, 0.16547804187493600915, 0.16767361861725019322,
165  0.16988540130252766513, 0.17211353531532005700, 0.17435816917135348788, 0.17661945459049489581,
166  0.17889754657247831241, 0.18119260347549629488, 0.18350478709776746150, 0.18583426276219711495,
167  0.18818119940425430485, 0.19054576966319540013, 0.19292814997677133873, 0.19532852067956322315,
168  0.19774706610509886464, 0.20018397469191127727, 0.20263943909370901930, 0.20511365629383770880,
169  0.20760682772422204205, 0.21011915938898825914, 0.21265086199297827522, 0.21520215107537867786,
170  0.21777324714870053264, 0.22036437584335949720, 0.22297576805812018050, 0.22560766011668406495,
171  0.22826029393071670664, 0.23093391716962742173, 0.23362878343743333945, 0.23634515245705964715,
172  0.23908329026244917002, 0.24184346939887722761, 0.24462596913189210901, 0.24743107566532763894,
173  0.25025908236886230967, 0.25311029001562948171, 0.25598500703041538015, 0.25888354974901621678,
174  0.26180624268936295243, 0.26475341883506220209, 0.26772541993204481808, 0.27072259679906003167,
175  0.27374530965280298302, 0.27679392844851734458, 0.27986883323697289920, 0.28297041453878076010,
176  0.28609907373707684673, 0.28925522348967773308, 0.29243928816189258772, 0.29565170428126120948,
177  0.29889292101558177099, 0.30216340067569352897, 0.30546361924459023541, 0.30879406693456016794,
178  0.31215524877417956945, 0.31554768522712893632, 0.31897191284495723773, 0.32242848495608914289,
179  0.32591797239355619822, 0.32944096426413633091, 0.33299806876180896713, 0.33658991402867758144,
180  0.34021714906678004560, 0.34388044470450243010, 0.34758049462163698567, 0.35131801643748334681,
181  0.35509375286678745925, 0.35890847294874976196, 0.36276297335481777335, 0.36665807978151414890,
182  0.37059464843514599421, 0.37457356761590215193, 0.37859575940958081092, 0.38266218149600982112,
183  0.38677382908413768115, 0.39093173698479710717, 0.39513698183329015336, 0.39939068447523107877,
184  0.40369401253053026739, 0.40804818315203238238, 0.41245446599716116772, 0.41691418643300289465,
185  0.42142872899761659635, 0.42599954114303435739, 0.43062813728845883923, 0.43531610321563659758,
186  0.44006510084235387501, 0.44487687341454851593, 0.44975325116275498919, 0.45469615747461548049,
187  0.45970761564213768669, 0.46478975625042618067, 0.46994482528395999841, 0.47517519303737738299,
188  0.48048336393045423016, 0.48587198734188493564, 0.49134386959403255500, 0.49690198724154955294,
189  0.50254950184134769289, 0.50828977641064283495, 0.51412639381474855788, 0.52006317736823356823,
190  0.52610421398361972602, 0.53225388026304326945, 0.53851687200286186590, 0.54489823767243963663,
191  0.55140341654064131685, 0.55803828226258748140, 0.56480919291240022434, 0.57172304866482579008,
192  0.57878735860284503057, 0.58601031847726802755, 0.59340090169173341521, 0.60096896636523224742,
193  0.60872538207962206507, 0.61668218091520762326, 0.62485273870366592605, 0.63325199421436607968,
194  0.64189671642726607018, 0.65080583341457104881, 0.66000084107899974178, 0.66950631673192477684,
195  0.67935057226476538741, 0.68956649611707798890, 0.70019265508278816709, 0.71127476080507597882,
196  0.72286765959357200702, 0.73503809243142351530, 0.74786862198519510742, 0.76146338884989624862,
197  0.77595685204011559675, 0.79152763697249565519, 0.80842165152300838005, 0.82699329664305033399,
198  0.84778550062398962096, 0.87170433238120363669, 0.90046992992574643800, 0.93814368086217467916,
199  1
200  };
201 
202  template <class RealType = double>
204  {
205  template <class Engine>
206  RealType operator()(Engine& eng)
207  {
208  const double* const table_x = exponential_table<double>::table_x;
209  const double* const table_y = exponential_table<double>::table_y;
210  RealType shift(0);
211  for (;;)
212  {
213  std::pair<RealType, int> vals = generate_int_float_pair_v165<RealType, 8>(eng);
214  int i = vals.second;
215  RealType x = vals.first * RealType(table_x[i]);
216  if (x < RealType(table_x[i + 1]))
217  return shift + x;
218  // For i=0 we need to generate from the tail, but because this is an exponential
219  // distribution, the tail looks exactly like the body, so we can simply repeat with a
220  // shift:
221  if (i == 0)
222  shift += RealType(table_x[1]);
223  else
224  {
225  RealType y01 = uniform_01<RealType>()(eng);
226  RealType y = RealType(table_y[i]) + y01 * RealType(table_y[i + 1] - table_y[i]);
227 
228  // All we care about is whether these are < or > 0; these values are equal to
229  // (lbound) or proportional to (ubound) `y` minus the lower/upper bound.
230  RealType y_above_ubound = RealType(table_x[i] - table_x[i + 1]) * y01 - (RealType(table_x[i]) - x),
231  y_above_lbound = y - (RealType(table_y[i + 1]) + (RealType(table_x[i + 1]) - x) * RealType(table_y[i + 1]));
232 
233  if (y_above_ubound < 0 // if above the upper bound reject immediately
234  && (y_above_lbound < 0 // If below the lower bound accept immediately
235  || y < f(x) // Otherwise it's between the bounds and we need a full check
236  ))
237  {
238  return x + shift;
239  }
240  }
241  }
242  }
243  static RealType f(RealType x)
244  {
245  using std::exp;
246  return exp(-x);
247  }
248  };
249 
250  } // namespace detail
251 
266  template <class RealType = double>
268  {
269  public:
270  typedef RealType input_type;
271  typedef RealType result_type;
272 
274  {
275  public:
277 
283  param_type(RealType lambda_arg = RealType(1.0))
284  : _lambda(lambda_arg) { BOOST_ASSERT(_lambda > RealType(0)); }
285 
287  RealType lambda() const { return _lambda; }
288 
291  {
292  os << parm._lambda;
293  return os;
294  }
295 
298  {
299  is >> parm._lambda;
300  return is;
301  }
302 
305  {
306  return lhs._lambda == rhs._lambda;
307  }
308 
310  BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
311 
312  private:
313  RealType _lambda;
314  };
315 
321  explicit exponential_distribution_v165(RealType lambda_arg = RealType(1.0))
322  : _lambda(lambda_arg) { BOOST_ASSERT(_lambda > RealType(0)); }
323 
328  : _lambda(parm.lambda()) {}
329 
330  // compiler-generated copy ctor and assignment operator are fine
331 
333  RealType lambda() const { return _lambda; }
334 
336  RealType min BOOST_PREVENT_MACRO_SUBSTITUTION() const
337  {
338  return RealType(0);
339  }
341  RealType max BOOST_PREVENT_MACRO_SUBSTITUTION() const
342  {
343  return (std::numeric_limits<RealType>::infinity)();
344  }
345 
347  param_type param() const { return param_type(_lambda); }
349  void param(const param_type& parm) { _lambda = parm.lambda(); }
350 
355  void reset() {}
356 
361  template <class Engine>
362  result_type operator()(Engine& eng) const
363  {
365  return impl(eng) / _lambda;
366  }
367 
372  template <class Engine>
373  result_type operator()(Engine& eng, const param_type& parm) const
374  {
375  return exponential_distribution_v165(parm)(eng);
376  }
377 
380  {
381  os << ed._lambda;
382  return os;
383  }
384 
387  {
388  is >> ed._lambda;
389  return is;
390  }
391 
397  {
398  return lhs._lambda == rhs._lambda;
399  }
400 
405  BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(exponential_distribution_v165)
406 
407  private:
408  result_type _lambda;
409  };
410 
411 } // namespace random
412 
413 } // namespace boost
414 
415 #endif // BOOST_165_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP
BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(exponential_distribution_v165, lhs, rhs)
BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, exponential_distribution_v165, ed)
BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, exponential_distribution_v165, ed)
result_type operator()(Engine &eng, const param_type &parm) const
exponential_distribution_v165(RealType lambda_arg=RealType(1.0))