/home/docs/checkouts/readthedocs.org/user_builds/advanced-micro-devices-rocrand/checkouts/docs-6.1.2/library/include/rocrand/rocrand_poisson.h Source File

/home/docs/checkouts/readthedocs.org/user_builds/advanced-micro-devices-rocrand/checkouts/docs-6.1.2/library/include/rocrand/rocrand_poisson.h Source File#

API library: /home/docs/checkouts/readthedocs.org/user_builds/advanced-micro-devices-rocrand/checkouts/docs-6.1.2/library/include/rocrand/rocrand_poisson.h Source File
API library
rocrand_poisson.h
1 // Copyright (c) 2017-2022 Advanced Micro Devices, Inc. All rights reserved.
2 //
3 // Permission is hereby granted, free of charge, to any person obtaining a copy
4 // of this software and associated documentation files (the "Software"), to deal
5 // in the Software without restriction, including without limitation the rights
6 // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
7 // copies of the Software, and to permit persons to whom the Software is
8 // furnished to do so, subject to the following conditions:
9 //
10 // The above copyright notice and this permission notice shall be included in
11 // all copies or substantial portions of the Software.
12 //
13 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
14 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
15 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
16 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
17 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
18 // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
19 // THE SOFTWARE.
20 
21 #ifndef ROCRAND_POISSON_H_
22 #define ROCRAND_POISSON_H_
23 
24 #ifndef FQUALIFIERS
25 #define FQUALIFIERS __forceinline__ __device__
26 #endif // FQUALIFIERS
27 
33 #include <math.h>
34 
35 #include "rocrand/rocrand_lfsr113.h"
36 #include "rocrand/rocrand_mrg31k3p.h"
37 #include "rocrand/rocrand_mrg32k3a.h"
38 #include "rocrand/rocrand_mtgp32.h"
39 #include "rocrand/rocrand_philox4x32_10.h"
40 #include "rocrand/rocrand_scrambled_sobol32.h"
41 #include "rocrand/rocrand_scrambled_sobol64.h"
42 #include "rocrand/rocrand_sobol32.h"
43 #include "rocrand/rocrand_sobol64.h"
44 #include "rocrand/rocrand_threefry2x32_20.h"
45 #include "rocrand/rocrand_threefry2x64_20.h"
46 #include "rocrand/rocrand_threefry4x32_20.h"
47 #include "rocrand/rocrand_threefry4x64_20.h"
48 #include "rocrand/rocrand_xorwow.h"
49 
50 #include "rocrand/rocrand_normal.h"
51 #include "rocrand/rocrand_uniform.h"
52 
53 namespace rocrand_device {
54 namespace detail {
55 
56 constexpr double lambda_threshold_small = 64.0;
57 constexpr double lambda_threshold_huge = 4000.0;
58 
59 template<class State, typename Result_Type = unsigned int>
60 FQUALIFIERS Result_Type poisson_distribution_small(State& state, double lambda)
61 {
62  // Knuth's method
63 
64  const double limit = exp(-lambda);
65  Result_Type k = 0;
66  double product = 1.0;
67 
68  do
69  {
70  k++;
71  product *= rocrand_uniform_double(state);
72  }
73  while (product > limit);
74 
75  return k - 1;
76 }
77 
79 double lgamma_approx(const double x)
80 {
81  // Lanczos approximation (g = 7, n = 9)
82 
83  const double z = x - 1.0;
84 
85  const int g = 7;
86  const int n = 9;
87  const double coefs[n] = {
88  0.99999999999980993227684700473478,
89  676.520368121885098567009190444019,
90  -1259.13921672240287047156078755283,
91  771.3234287776530788486528258894,
92  -176.61502916214059906584551354,
93  12.507343278686904814458936853,
94  -0.13857109526572011689554707,
95  9.984369578019570859563e-6,
96  1.50563273514931155834e-7
97  };
98  double sum = 0.0;
99  #pragma unroll
100  for (int i = n - 1; i > 0; i--)
101  {
102  sum += coefs[i] / (z + i);
103  }
104  sum += coefs[0];
105 
106  const double log_sqrt_2_pi = 0.9189385332046727418;
107  const double e = 2.718281828459045090796;
108  return (log_sqrt_2_pi + log(sum) - g) + (z + 0.5) * log((z + g + 0.5) / e);
109 }
110 
111 template<class State, typename Result_Type = unsigned int>
112 FQUALIFIERS Result_Type poisson_distribution_large(State& state, double lambda)
113 {
114  // Rejection method PA, A. C. Atkinson
115 
116  const double c = 0.767 - 3.36 / lambda;
117  const double beta = ROCRAND_PI_DOUBLE / sqrt(3.0 * lambda);
118  const double alpha = beta * lambda;
119  const double k = log(c) - lambda - log(beta);
120  const double log_lambda = log(lambda);
121 
122  while (true)
123  {
124  const double u = rocrand_uniform_double(state);
125  const double x = (alpha - log((1.0 - u) / u)) / beta;
126  const double n = floor(x + 0.5);
127  if (n < 0)
128  {
129  continue;
130  }
131  const double v = rocrand_uniform_double(state);
132  const double y = alpha - beta * x;
133  const double t = 1.0 + exp(y);
134  const double lhs = y + log(v / (t * t));
135  const double rhs = k + n * log_lambda - lgamma_approx(n + 1.0);
136  if (lhs <= rhs)
137  {
138  return static_cast<Result_Type>(n);
139  }
140  }
141 }
142 
143 template<class State, typename Result_Type = unsigned int>
144 FQUALIFIERS Result_Type poisson_distribution_huge(State& state, double lambda)
145 {
146  // Approximate Poisson distribution with normal distribution
147 
148  const double n = rocrand_normal_double(state);
149  return static_cast<Result_Type>(round(sqrt(lambda) * n + lambda));
150 }
151 
152 template<class State, typename Result_Type = unsigned int>
153 FQUALIFIERS Result_Type poisson_distribution(State& state, double lambda)
154 {
155  if (lambda < lambda_threshold_small)
156  {
157  return poisson_distribution_small<State, Result_Type>(state, lambda);
158  }
159  else if (lambda <= lambda_threshold_huge)
160  {
161  return poisson_distribution_large<State, Result_Type>(state, lambda);
162  }
163  else
164  {
165  return poisson_distribution_huge<State, Result_Type>(state, lambda);
166  }
167 }
168 
169 template<class State, typename Result_Type = unsigned int>
170 FQUALIFIERS Result_Type poisson_distribution_itr(State& state, double lambda)
171 {
172  // Algorithm ITR
173  // George S. Fishman
174  // Discrete-Event Simulation: Modeling, Programming, and Analysis
175  // p. 333
176  double L;
177  double x = 1.0;
178  double y = 1.0;
179  Result_Type k = 0;
180  int pow = 0;
181  // Algorithm ITR uses u from (0, 1) and uniform_double returns (0, 1]
182  // Change u to ensure that 1 is never generated,
183  // otherwise the inner loop never ends.
184  double u = rocrand_uniform_double(state) - ROCRAND_2POW32_INV_DOUBLE / 2.0;
185  double upow = pow + 500.0;
186  double ex = exp(-500.0);
187  do{
188  if (lambda > upow)
189  L = ex;
190  else
191  L = exp((double)(pow - lambda));
192 
193  x *= L;
194  y *= L;
195  pow += 500;
196  while (u > y)
197  {
198  k++;
199  x *= ((double)lambda / (double) k);
200  y += x;
201  }
202  } while((double)pow < lambda);
203  return k;
204 }
205 
206 template<class State, typename Result_Type = unsigned int>
207 FQUALIFIERS Result_Type poisson_distribution_inv(State& state, double lambda)
208 {
209  if (lambda < 1000.0)
210  {
211  return poisson_distribution_itr<State, Result_Type>(state, lambda);
212  }
213  else
214  {
215  return poisson_distribution_huge<State, Result_Type>(state, lambda);
216  }
217 }
218 
219 } // end namespace detail
220 } // end namespace rocrand_device
221 
233 #ifndef ROCRAND_DETAIL_PHILOX_BM_NOT_IN_STATE
235 unsigned int rocrand_poisson(rocrand_state_philox4x32_10 * state, double lambda)
236 {
237  return rocrand_device::detail::poisson_distribution<rocrand_state_philox4x32_10*, unsigned int>(
238  state,
239  lambda);
240 }
241 
254 uint4 rocrand_poisson4(rocrand_state_philox4x32_10 * state, double lambda)
255 {
256  return uint4{
257  rocrand_device::detail::poisson_distribution<rocrand_state_philox4x32_10*, unsigned int>(
258  state,
259  lambda),
260  rocrand_device::detail::poisson_distribution<rocrand_state_philox4x32_10*, unsigned int>(
261  state,
262  lambda),
263  rocrand_device::detail::poisson_distribution<rocrand_state_philox4x32_10*, unsigned int>(
264  state,
265  lambda),
266  rocrand_device::detail::poisson_distribution<rocrand_state_philox4x32_10*, unsigned int>(
267  state,
268  lambda)};
269 }
270 #endif // ROCRAND_DETAIL_PHILOX_BM_NOT_IN_STATE
271 
283 #ifndef ROCRAND_DETAIL_MRG31K3P_BM_NOT_IN_STATE
284 FQUALIFIERS unsigned int rocrand_poisson(rocrand_state_mrg31k3p* state, double lambda)
285 {
286  return rocrand_device::detail::poisson_distribution<rocrand_state_mrg31k3p*, unsigned int>(
287  state,
288  lambda);
289 }
290 #endif // ROCRAND_DETAIL_MRG31K3P_BM_NOT_IN_STATE
291 
303 #ifndef ROCRAND_DETAIL_MRG32K3A_BM_NOT_IN_STATE
305 unsigned int rocrand_poisson(rocrand_state_mrg32k3a * state, double lambda)
306 {
307  return rocrand_device::detail::poisson_distribution<rocrand_state_mrg32k3a*, unsigned int>(
308  state,
309  lambda);
310 }
311 #endif // ROCRAND_DETAIL_MRG32K3A_BM_NOT_IN_STATE
312 
324 #ifndef ROCRAND_DETAIL_XORWOW_BM_NOT_IN_STATE
326 unsigned int rocrand_poisson(rocrand_state_xorwow * state, double lambda)
327 {
328  return rocrand_device::detail::poisson_distribution<rocrand_state_xorwow*, unsigned int>(
329  state,
330  lambda);
331 }
332 #endif // ROCRAND_DETAIL_XORWOW_BM_NOT_IN_STATE
333 
346 unsigned int rocrand_poisson(rocrand_state_mtgp32 * state, double lambda)
347 {
348  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_mtgp32*, unsigned int>(
349  state,
350  lambda);
351 }
352 
365 unsigned int rocrand_poisson(rocrand_state_sobol32 * state, double lambda)
366 {
367  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_sobol32*, unsigned int>(
368  state,
369  lambda);
370 }
371 
384 unsigned int rocrand_poisson(rocrand_state_scrambled_sobol32* state, double lambda)
385 {
386  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_scrambled_sobol32*,
387  unsigned int>(state, lambda);
388 }
389 
402 unsigned long long int rocrand_poisson(rocrand_state_sobol64* state, double lambda)
403 {
404  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_sobol64*,
405  unsigned long long int>(state, lambda);
406 }
407 
420 unsigned long long int rocrand_poisson(rocrand_state_scrambled_sobol64* state, double lambda)
421 {
422  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_scrambled_sobol64*,
423  unsigned long long int>(state, lambda);
424 }
425 
438 unsigned int rocrand_poisson(rocrand_state_lfsr113* state, double lambda)
439 {
440  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_lfsr113*, unsigned int>(
441  state,
442  lambda);
443 }
444 
456 FQUALIFIERS unsigned int rocrand_poisson(rocrand_state_threefry2x32_20* state, double lambda)
457 {
458  return rocrand_device::detail::poisson_distribution_inv(state, lambda);
459 }
460 
472 FQUALIFIERS unsigned int rocrand_poisson(rocrand_state_threefry2x64_20* state, double lambda)
473 {
474  return rocrand_device::detail::poisson_distribution_inv(state, lambda);
475 }
476 
488 FQUALIFIERS unsigned int rocrand_poisson(rocrand_state_threefry4x32_20* state, double lambda)
489 {
490  return rocrand_device::detail::poisson_distribution_inv(state, lambda);
491 }
492 
504 FQUALIFIERS unsigned int rocrand_poisson(rocrand_state_threefry4x64_20* state, double lambda)
505 {
506  return rocrand_device::detail::poisson_distribution_inv(state, lambda);
507 }
508  // end of group rocranddevice
510 
511 #endif // ROCRAND_POISSON_H_
FQUALIFIERS double rocrand_normal_double(rocrand_state_philox4x32_10 *state)
Returns a normally distributed double value.
Definition: rocrand_normal.h:456
FQUALIFIERS double rocrand_uniform_double(rocrand_state_philox4x32_10 *state)
Returns a uniformly distributed random double value from (0; 1] range.
Definition: rocrand_uniform.h:303
FQUALIFIERS uint4 rocrand_poisson4(rocrand_state_philox4x32_10 *state, double lambda)
Returns four Poisson-distributed unsigned int values using Philox generator.
Definition: rocrand_poisson.h:254
FQUALIFIERS unsigned int rocrand_poisson(rocrand_state_philox4x32_10 *state, double lambda)
Returns a Poisson-distributed unsigned int using Philox generator.
Definition: rocrand_poisson.h:235
#define FQUALIFIERS
Shorthand for commonly used function qualifiers.
Definition: rocrand_uniform.h:31