HDK
UT_Poisson.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: UT_Poisson.h (UT Library, C++)
7  *
8  * COMMENTS: Generating discrete random Poisson variates.
9  */
10
11 #ifndef __UT_Poisson__
12 #define __UT_Poisson__
13
14 #include <SYS/SYS_Floor.h>
15 #include <SYS/SYS_Math.h>
16 #include <SYS/SYS_Types.h>
17
19 {
20 public:
21  template<typename T>
22  static int64 GetPoisson(fpreal mean, T &generator, bool limit_range=false, int limit_min=-1, int limit_max=-1)
23  {
24  if (mean <= 0.0f)
25  return 0;
26
27  if (limit_range)
28  {
29  if (limit_min > limit_max)
30  {
31  int tmp = limit_min;
32  limit_min = limit_max;
33  limit_max = tmp;
34  }
35
36  if (limit_min < 0)
37  limit_min = 0;
38  if (limit_max < 0)
39  return 0;
40  }
41
42  int max_tries = 1000;
43
44  // For small mean, just sum the probabilities until they're
45  // greater than or equal to a uniform random number.
46  const float smallmeanthreshold = 19;
47  // (50-19)/sqrt(19) = 7.1 standard deviations, so we'll never hit it
48  // unless things have gone wrong or maybe if fpreal becomes 32-bit.
49  const int64 smallmeanmaxi = 50;
50  if (mean < smallmeanthreshold)
51  {
52  int tries = 0;
53  while (tries < max_tries)
54  {
55  tries++;
56
57  const float u = generator.getFloat();
58  // Exact equality check, just to avoid ln(0)
59  if (u == 0)
60  {
61  if (limit_range && (0 < limit_min || 0 > limit_max))
62  continue;
63  else
64  return 0;
65  }
66
67  fpreal r0 = u * SYSexp(mean);
68  if (r0 <= 1.0)
69  {
70  if (limit_range && (0 < limit_min || 0 > limit_max))
71  continue;
72  else
73  return 0;
74  }
75
76  fpreal f = 1.0 + mean;
77  if (r0 <= f)
78  {
79  if (limit_range && (1 < limit_min || 1 > limit_max))
80  continue;
81  else
82  return 1;
83  }
84
85  fpreal term = 0.5*mean*mean;
86  f += term;
87  int64 i = 2;
88  while (r0 > f && i < smallmeanmaxi && (!limit_range || i <= limit_max))
89  {
90  ++i;
91  term *= (mean / i);
92  f += term;
93  }
94
95  if (limit_range && (i < limit_min || i > limit_max))
96  continue;
97
98  return i;
99  }
100
101  return limit_min; // or max?
102  }
103
104  // TODO: Which SYSfloor to use?
105  // fpreal or float? why?
106  const int64 mode = int64(mean);
107  const fpreal a = mean + 0.5;
108  const int64 k = int64(SYSfloor(a - SYSsqrt(2.0*a))) + 1;
109  fpreal f_k = 1;
110  if (k == mode)
111  {
112  // do nothing. It's one.
113  }
114  else if (k > mode)
115  {
116  int64 curr = k;
117  while (curr > mode)
118  {
119  f_k *= (mean / fpreal(curr));
120  curr--;
121  }
122  }
123  else // (k < mode)
124  {
125  int64 curr = mode;
126  while (curr > k)
127  {
128  f_k *= (fpreal(curr) / mean);
129  curr--;
130  }
131  }
132
133  const fpreal s = (a - k)*SYSsqrt(f_k);
134  // Now we have (s) and (a)
135  int64 result = 0;
136  fpreal x, y, f_x;
137  int tries = 0;
138  while (!limit_range || tries < max_tries)
139  {
140  tries++;
141
142  // Generate u and v
143  const float u = generator.getFloat();
144  const float v = generator.getFloat();
145  // set X
146  x = a + (s * (2 * v - 1) / u);
147  // reject all negative results
148  if (x < 0)
149  continue;
150  // set K (= result)
151  result = int64(SYSfloor(x));
152
153  // if there is a range limit, and it falls outside of the range
154  if (limit_range && (result < limit_min || result > limit_max))
155  continue;
156
157  // set Y
158  y = u * u;
159  // evaluate f_x and compare against y
160  f_x = 1;
161  if (result == mode)
162  {
163  return result;
164  }
165  else if (result > mode)
166  {
167  int64 curr = result;
168  bool stopped = false;
169  while (curr > mode)
170  {
171  f_x *= (mean / fpreal(curr));
172  if (y > f_x)
173  {
174  stopped = true;
175  break;
176  }
177  curr--;
178  }
179  if (stopped)
180  continue;
181  }
182  else // (result < mode)
183  {
184  int64 curr = mode;
185  bool stopped = false;
186  while (curr > result)
187  {
188  f_x *= (fpreal(curr) / mean);
189  if (y > f_x)
190  {
191  stopped = true;
192  break;
193  }
194  curr--;
195  }
196  if (stopped)
197  continue;
198  }
199
200  if (y <= f_x)
201  return result;
202  }
203
204  return limit_min; // or max?
205  }
206
207 };
208
209 #endif
static int64 GetPoisson(fpreal mean, T &generator, bool limit_range=false, int limit_min=-1, int limit_max=-1)
Definition: UT_Poisson.h:22
GLint GLenum GLint x
Definition: glcorearb.h:409
GLuint64EXT * result
Definition: glew.h:14311
const GLdouble * v
Definition: glcorearb.h:837
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
long long int64
Definition: SYS_Types.h:116
SYS_API fpreal32 SYSfloor(fpreal32 val)
GLenum mode
Definition: glcorearb.h:99
fpreal64 fpreal
Definition: SYS_Types.h:277
GLfloat f
Definition: glcorearb.h:1926
GLdouble s
Definition: glew.h:1395
GLint y
Definition: glcorearb.h:103