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

API library

• library
• include
• rocrand

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 

78 FQUALIFIERS

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

234 FQUALIFIERS

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 

253 FQUALIFIERS

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

304 FQUALIFIERS

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

325 FQUALIFIERS

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 

345 FQUALIFIERS

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 

364 FQUALIFIERS

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 

383 FQUALIFIERS

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 

401 FQUALIFIERS

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 

419 FQUALIFIERS

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 

437 FQUALIFIERS

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_

rocrand_normal_double

FQUALIFIERS double rocrand_normal_double(rocrand_state_philox4x32_10 *state)

Returns a normally distributed double value.

Definition: rocrand_normal.h:456

rocrand_uniform_double

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

rocrand_poisson4

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

rocrand_poisson

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

FQUALIFIERS

#define FQUALIFIERS

Shorthand for commonly used function qualifiers.

Definition: rocrand_uniform.h:31

---

Generated by  1.9.1