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

• library
• include
• rocrand

rocrand_poisson.h

1 // Copyright (c) 2017-2024 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 

29 #include <math.h>

30 

31 #include "rocrand/rocrand_lfsr113.h"

32 #include "rocrand/rocrand_mrg31k3p.h"

33 #include "rocrand/rocrand_mrg32k3a.h"

34 #include "rocrand/rocrand_mtgp32.h"

35 #include "rocrand/rocrand_philox4x32_10.h"

36 #include "rocrand/rocrand_scrambled_sobol32.h"

37 #include "rocrand/rocrand_scrambled_sobol64.h"

38 #include "rocrand/rocrand_sobol32.h"

39 #include "rocrand/rocrand_sobol64.h"

40 #include "rocrand/rocrand_threefry2x32_20.h"

41 #include "rocrand/rocrand_threefry2x64_20.h"

42 #include "rocrand/rocrand_threefry4x32_20.h"

43 #include "rocrand/rocrand_threefry4x64_20.h"

44 #include "rocrand/rocrand_xorwow.h"

45 

46 #include "rocrand/rocrand_normal.h"

47 #include "rocrand/rocrand_uniform.h"

48 

49 namespace rocrand_device {

50 namespace detail {

51 

52 constexpr double lambda_threshold_small = 64.0;

53 constexpr double lambda_threshold_huge = 4000.0;

54 

55 template<class State, typename Result_Type = unsigned int>

56 __forceinline__ __device__ __host__ Result_Type poisson_distribution_small(State& state,

57  double lambda)

58 {

59  // Knuth's method

60 

61  const double limit = exp(-lambda);

62  Result_Type k = 0;

63  double product = 1.0;

64 

65  do

66  {

67  k++;

68  product *= rocrand_uniform_double(state);

69  }

70  while (product > limit);

71 

72  return k - 1;

73 }

74 

75 __forceinline__ __device__ __host__ double lgamma_approx(const double x)

76 {

77  // Lanczos approximation (g = 7, n = 9)

78 

79  const double z = x - 1.0;

80 

81  const int g = 7;

82  const int n = 9;

83  const double coefs[n] = {

84  0.99999999999980993227684700473478,

85  676.520368121885098567009190444019,

86  -1259.13921672240287047156078755283,

87  771.3234287776530788486528258894,

88  -176.61502916214059906584551354,

89  12.507343278686904814458936853,

90  -0.13857109526572011689554707,

91  9.984369578019570859563e-6,

92  1.50563273514931155834e-7

93  };

94  double sum = 0.0;

95  #pragma unroll

96  for (int i = n - 1; i > 0; i--)

97  {

98  sum += coefs[i] / (z + i);

99  }

100  sum += coefs[0];

101 

102  const double log_sqrt_2_pi = 0.9189385332046727418;

103  const double e = 2.718281828459045090796;

104  return (log_sqrt_2_pi + log(sum) - g) + (z + 0.5) * log((z + g + 0.5) / e);

105 }

106 

107 template<class State, typename Result_Type = unsigned int>

108 __forceinline__ __device__ __host__ Result_Type poisson_distribution_large(State& state,

109  double lambda)

110 {

111  // Rejection method PA, A. C. Atkinson

112 

113  const double c = 0.767 - 3.36 / lambda;

114  const double beta = ROCRAND_PI_DOUBLE / sqrt(3.0 * lambda);

115  const double alpha = beta * lambda;

116  const double k = log(c) - lambda - log(beta);

117  const double log_lambda = log(lambda);

118 

119  while (true)

120  {

121  const double u = rocrand_uniform_double(state);

122  const double x = (alpha - log((1.0 - u) / u)) / beta;

123  const double n = floor(x + 0.5);

124  if (n < 0)

125  {

126  continue;

127  }

128  const double v = rocrand_uniform_double(state);

129  const double y = alpha - beta * x;

130  const double t = 1.0 + exp(y);

131  const double lhs = y + log(v / (t * t));

132  const double rhs = k + n * log_lambda - lgamma_approx(n + 1.0);

133  if (lhs <= rhs)

134  {

135  return static_cast<Result_Type>(n);

136  }

137  }

138 }

139 

140 template<class State, typename Result_Type = unsigned int>

141 __forceinline__ __device__ __host__ Result_Type poisson_distribution_huge(State& state,

142  double lambda)

143 {

144  // Approximate Poisson distribution with normal distribution

145 

146  const double n = rocrand_normal_double(state);

147  return static_cast<Result_Type>(round(sqrt(lambda) * n + lambda));

148 }

149 

150 template<class State, typename Result_Type = unsigned int>

151 __forceinline__ __device__ __host__ Result_Type poisson_distribution(State& state, double lambda)

152 {

153  if (lambda < lambda_threshold_small)

154  {

155  return poisson_distribution_small<State, Result_Type>(state, lambda);

156  }

157  else if (lambda <= lambda_threshold_huge)

158  {

159  return poisson_distribution_large<State, Result_Type>(state, lambda);

160  }

161  else

162  {

163  return poisson_distribution_huge<State, Result_Type>(state, lambda);

164  }

165 }

166 

167 template<class State, typename Result_Type = unsigned int>

168 __forceinline__ __device__ __host__ Result_Type poisson_distribution_itr(State& state,

169  double lambda)

170 {

171  // Algorithm ITR

172  // George S. Fishman

173  // Discrete-Event Simulation: Modeling, Programming, and Analysis

174  // p. 333

175  double L;

176  double x = 1.0;

177  double y = 1.0;

178  Result_Type k = 0;

179  int pow = 0;

180  // Algorithm ITR uses u from (0, 1) and uniform_double returns (0, 1]

181  // Change u to ensure that 1 is never generated,

182  // otherwise the inner loop never ends.

183  double u = rocrand_uniform_double(state) - ROCRAND_2POW32_INV_DOUBLE / 2.0;

184  double upow = pow + 500.0;

185  double ex = exp(-500.0);

186  do{

187  if (lambda > upow)

188  L = ex;

189  else

190  L = exp((double)(pow - lambda));

191 

192  x *= L;

193  y *= L;

194  pow += 500;

195  while (u > y)

196  {

197  k++;

198  x *= ((double)lambda / (double) k);

199  y += x;

200  }

201  } while((double)pow < lambda);

202  return k;

203 }

204 

205 template<class State, typename Result_Type = unsigned int>

206 __forceinline__ __device__ __host__ Result_Type poisson_distribution_inv(State& state,

207  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_BM_NOT_IN_STATE

234 __forceinline__ __device__ __host__ unsigned int rocrand_poisson(rocrand_state_philox4x32_10* state,

235  double lambda)

236 {

237  return rocrand_device::detail::poisson_distribution<rocrand_state_philox4x32_10*, unsigned int>(

238  state,

239  lambda);

240 }

241 

253 __forceinline__ __device__ __host__ uint4 rocrand_poisson4(rocrand_state_philox4x32_10* state,

254  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_BM_NOT_IN_STATE

271 

283 #ifndef ROCRAND_DETAIL_BM_NOT_IN_STATE

284 __forceinline__ __device__ __host__ unsigned int rocrand_poisson(rocrand_state_mrg31k3p* state,

285  double lambda)

286 {

287  return rocrand_device::detail::poisson_distribution<rocrand_state_mrg31k3p*, unsigned int>(

288  state,

289  lambda);

290 }

291 #endif // ROCRAND_DETAIL_BM_NOT_IN_STATE

292 

304 #ifndef ROCRAND_DETAIL_BM_NOT_IN_STATE

305 __forceinline__ __device__ __host__ unsigned int rocrand_poisson(rocrand_state_mrg32k3a* state,

306  double lambda)

307 {

308  return rocrand_device::detail::poisson_distribution<rocrand_state_mrg32k3a*, unsigned int>(

309  state,

310  lambda);

311 }

312 #endif // ROCRAND_DETAIL_BM_NOT_IN_STATE

313 

325 #ifndef ROCRAND_DETAIL_BM_NOT_IN_STATE

326 __forceinline__ __device__ __host__ unsigned int rocrand_poisson(rocrand_state_xorwow* state,

327  double lambda)

328 {

329  return rocrand_device::detail::poisson_distribution<rocrand_state_xorwow*, unsigned int>(

330  state,

331  lambda);

332 }

333 #endif // ROCRAND_DETAIL_BM_NOT_IN_STATE

334 

346 __forceinline__ __device__ __host__ unsigned int rocrand_poisson(rocrand_state_mtgp32* state,

347  double lambda)

348 {

349  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_mtgp32*, unsigned int>(

350  state,

351  lambda);

352 }

353 

365 __forceinline__ __device__ __host__ unsigned int rocrand_poisson(rocrand_state_sobol32* state,

366  double lambda)

367 {

368  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_sobol32*, unsigned int>(

369  state,

370  lambda);

371 }

372 

384 __forceinline__ __device__ __host__ unsigned int

385  rocrand_poisson(rocrand_state_scrambled_sobol32* state, double lambda)

386 {

387  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_scrambled_sobol32*,

388  unsigned int>(state, lambda);

389 }

390 

402 __forceinline__ __device__ __host__ unsigned long long int

403  rocrand_poisson(rocrand_state_sobol64* state, double lambda)

404 {

405  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_sobol64*,

406  unsigned long long int>(state, lambda);

407 }

408 

420 __forceinline__ __device__ __host__ unsigned long long int

421  rocrand_poisson(rocrand_state_scrambled_sobol64* state, double lambda)

422 {

423  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_scrambled_sobol64*,

424  unsigned long long int>(state, lambda);

425 }

426 

438 __forceinline__ __device__ __host__ unsigned int rocrand_poisson(rocrand_state_lfsr113* state,

439  double lambda)

440 {

441  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_lfsr113*, unsigned int>(

442  state,

443  lambda);

444 }

445 

457 __forceinline__ __device__ __host__ unsigned int

458  rocrand_poisson(rocrand_state_threefry2x32_20* state, double lambda)

459 {

460  return rocrand_device::detail::poisson_distribution_inv(state, lambda);

461 }

462 

474 __forceinline__ __device__ __host__ unsigned int

475  rocrand_poisson(rocrand_state_threefry2x64_20* state, double lambda)

476 {

477  return rocrand_device::detail::poisson_distribution_inv(state, lambda);

478 }

479 

491 __forceinline__ __device__ __host__ unsigned int

492  rocrand_poisson(rocrand_state_threefry4x32_20* state, double lambda)

493 {

494  return rocrand_device::detail::poisson_distribution_inv(state, lambda);

495 }

496 

508 __forceinline__ __device__ __host__ unsigned int

509  rocrand_poisson(rocrand_state_threefry4x64_20* state, double lambda)

510 {

511  return rocrand_device::detail::poisson_distribution_inv(state, lambda);

512 }

513  // end of group rocranddevice

515 

516 #endif // ROCRAND_POISSON_H_

rocrand_poisson

__forceinline__ __device__ __host__ unsigned int rocrand_poisson(rocrand_state_philox4x32_10 *state, double lambda)

Returns a Poisson-distributed unsigned int using Philox generator.

Definition: rocrand_poisson.h:234

rocrand_poisson4

__forceinline__ __device__ __host__ uint4 rocrand_poisson4(rocrand_state_philox4x32_10 *state, double lambda)

Returns four Poisson-distributed unsigned int values using Philox generator.

Definition: rocrand_poisson.h:253

rocrand_uniform_double

__forceinline__ __device__ __host__ double rocrand_uniform_double(rocrand_state_philox4x32_10 *state)

Returns a uniformly distributed random double value from (0; 1] range.

Definition: rocrand_uniform.h:294

rocrand_normal_double

__forceinline__ __device__ __host__ double rocrand_normal_double(rocrand_state_philox4x32_10 *state)

Returns a normally distributed double value.

Definition: rocrand_normal.h:438