libStatGen Software  1
Random.cpp
1 /*
2  * Copyright (C) 2010 Regents of the University of Michigan
3  *
4  * This program is free software: you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation, either version 3 of the License, or
7  * (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program. If not, see <http://www.gnu.org/licenses/>.
16  */
17 
18 
19 //////////////////////////////////////////////////////////////////////////////
20 // This file includes code derived from the original Mersenne Twister Code
21 // by Makoto Matsumoto and Takuji Nishimura
22 // and is subject to their original copyright notice copied below:
23 //////////////////////////////////////////////////////////////////////////////
24 
25 //////////////////////////////////////////////////////////////////////////////
26 // COPYRIGHT NOTICE FOR MERSENNE TWISTER CODE
27 // Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
28 // All rights reserved.
29 //
30 // Redistribution and use in source and binary forms, with or without
31 // modification, are permitted provided that the following conditions
32 // are met:
33 //
34 // 1. Redistributions of source code must retain the above copyright
35 // notice, this list of conditions and the following disclaimer.
36 //
37 // 2. Redistributions in binary form must reproduce the above copyright
38 // notice, this list of conditions and the following disclaimer in the
39 // documentation and/or other materials provided with the distribution.
40 //
41 // 3. The names of its contributors may not be used to endorse or promote
42 // products derived from this software without specific prior written
43 // permission.
44 //
45 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
46 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
47 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
48 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
49 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
50 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
51 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
52 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
53 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
54 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
55 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
56 //
57 ///////////////////////////////////////////////////////////////////////////////
58 
59 #include "Random.h"
60 #include "MathConstant.h"
61 #include "Error.h"
62 
63 #include <math.h>
64 
65 //Constants used internally by Mersenne random number generator
66 #define MERSENNE_N 624
67 #define MERSENNE_M 397
68 
69 // constant vector a
70 #define MATRIX_A 0x9908b0dfUL
71 
72 // most significant w-r bits
73 #define UPPER_MASK 0x80000000UL
74 
75 // least significant r bits
76 #define LOWER_MASK 0x7fffffffUL
77 
78 
79 // Constants used internally by Park-Miller random generator
80 #define IA 16807
81 #define IM 2147483647
82 #define AM (1.0 / IM)
83 #define IQ 127773
84 #define IR 2836
85 #define NTAB 32
86 #define NDIV (1+(IM-1)/NTAB)
87 #define RNMX (1.0-EPS)
88 
89 Random::Random(long s)
90 {
91 #ifndef __NO_MERSENNE
92  mt = new unsigned long [MERSENNE_N];
93  mti = MERSENNE_N + 1;
94  mersenneMult = 1.0/4294967296.0;
95 #else
96  shuffler = new long [NTAB];
97 #endif
98  Reset(s);
99 }
100 
101 Random::~Random()
102 {
103 #ifndef __NO_MERSENNE
104  delete [] mt;
105 #else
106  delete [] shuffler;
107 #endif
108 }
109 
110 void Random::Reset(long s)
111 {
112  normSaved = 0;
113 
114 #ifndef __NO_MERSENNE
115  InitMersenne(s);
116 #else
117  // 'Continuous' Random Generator
118  if ((seed = s) < 1)
119  seed = s == 0 ? 1 : -s; // seed == 0 would be disastrous
120 
121  for (int j=NTAB+7; j>=0; j--) // Warm up and set shuffle table
122  {
123  long k = seed / IQ;
124  seed = IA * (seed - k * IQ) - IR * k;
125  if (seed < 0) seed += IM;
126  if (j < NTAB) shuffler[j] = seed;
127  }
128  last=shuffler[0];
129 #endif
130 }
131 
132 // initializes mt[MERSENNE_N] with a seed
133 void Random::InitMersenne(unsigned long s)
134 {
135  mt[0]= s & 0xffffffffUL;
136  for (mti = 1; mti < MERSENNE_N; mti++)
137  {
138  mt[mti] = (1812433253UL * (mt[mti-1] ^(mt[mti-1] >> 30)) + mti);
139  /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
140  /* In the previous versions, MSBs of the seed affect */
141  /* only MSBs of the array mt[]. */
142  /* 2002/01/09 modified by Makoto Matsumoto */
143 
144  mt[mti] &= 0xffffffffUL;
145  }
146 }
147 
148 int Random::Binary()
149 {
150  return Next() > 0.5 ? 1 : 0;
151 }
152 
153 #ifndef __NO_MERSENNE
154 
155 double Random::Next()
156 {
157  unsigned long y;
158 
159  // mag01[x] = x * MATRIX_A for x=0,1
160  static unsigned long mag01[2]={0x0UL, MATRIX_A};
161 
162  if (mti >= MERSENNE_N)
163  {
164  /* generate MERSENNE_N words at one time */
165  int kk;
166 
167  // If InitMersenne() has not been called, a default initial seed is used
168  if (mti == MERSENNE_N+1)
169  InitMersenne(5489UL);
170 
171  for (kk=0; kk < MERSENNE_N-MERSENNE_M; kk++)
172  {
173  y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
174  mt[kk] = mt[kk+MERSENNE_M] ^(y >> 1) ^ mag01[y & 0x1UL];
175  }
176 
177  for (; kk < MERSENNE_N-1; kk++)
178  {
179  y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
180  mt[kk] = mt[kk+(MERSENNE_M - MERSENNE_N)] ^(y >> 1) ^ mag01[y & 0x1UL];
181  }
182 
183  y = (mt[MERSENNE_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
184  mt[MERSENNE_N-1] = mt[MERSENNE_M-1] ^(y >> 1) ^ mag01[y & 0x1UL];
185 
186  mti = 0;
187  }
188 
189  y = mt[mti++];
190 
191  // Tempering
192  y ^= (y >> 11);
193  y ^= (y << 7) & 0x9d2c5680UL;
194  y ^= (y << 15) & 0xefc60000UL;
195  y ^= (y >> 18);
196 
197  return (mersenneMult *((double) y + 0.5));
198 }
199 
200 // Generates a random number on [0,0xffffffff]-interval
201 
202 unsigned long Random::NextInt()
203 {
204  unsigned long y;
205 
206  // mag01[x] = x * MATRIX_A for x=0,1
207  static unsigned long mag01[2]={0x0UL, MATRIX_A};
208 
209  if (mti >= MERSENNE_N)
210  {
211  /* generate MERSENNE_N words at one time */
212  int kk;
213 
214  // If InitMersenne() has not been called, a default initial seed is used
215  if (mti == MERSENNE_N + 1)
216  InitMersenne(5489UL);
217 
218  for (kk= 0; kk < MERSENNE_N - MERSENNE_M; kk++)
219  {
220  y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
221  mt[kk] = mt[kk+MERSENNE_M] ^(y >> 1) ^ mag01[y & 0x1UL];
222  }
223 
224  for (; kk< MERSENNE_N-1; kk++)
225  {
226  y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
227  mt[kk] = mt[kk+(MERSENNE_M - MERSENNE_N)] ^(y >> 1) ^ mag01[y & 0x1UL];
228  }
229 
230  y = (mt[MERSENNE_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
231  mt[MERSENNE_N-1] = mt[MERSENNE_M-1] ^(y >> 1) ^ mag01[y & 0x1UL];
232 
233  mti = 0;
234  }
235 
236  y = mt[mti++];
237 
238  // Tempering
239  y ^= (y >> 11);
240  y ^= (y << 7) & 0x9d2c5680UL;
241  y ^= (y << 15) & 0xefc60000UL;
242  y ^= (y >> 18);
243 
244  return y;
245 }
246 
247 #else
248 
249 double Random::Next()
250 {
251  // Compute seed = (IA * seed) % IM without overflows
252  // by Schrage's method
253  long k = seed / IQ;
254  seed = IA * (seed - k * IQ) - IR * k;
255  if (seed < 0) seed += IM;
256 
257  // Map to 0..NTAB-1
258  int j = last/NDIV;
259 
260  // Output value is shuffler[j], which is in turn replaced by seed
261  last = shuffler[j];
262  shuffler[j] = seed;
263 
264  // Map to 0.0 .. 1.0 excluding endpoints
265  double temp = AM * last;
266  if (temp > RNMX) return RNMX;
267  return temp;
268 }
269 
270 unsigned long Random::NextInt()
271 {
272  // Compute seed = (IA * seed) % IM without overflows
273  // by Schrage's method
274  long k = seed / IQ;
275  seed = IA * (seed - k * IQ) - IR * k;
276  if (seed < 0) seed += IM;
277 
278  // Map to 0..NTAB-1
279  int j = last/NDIV;
280 
281  // Output value is shuffler[j], which is in turn replaced by seed
282  last = shuffler[j];
283  shuffler[j] = seed;
284 
285  return last;
286 }
287 
288 #endif
289 
290 double Random::Normal()
291 {
292  double v1, v2, fac, rsq;
293 
294  if (!normSaved) // Do we need new numbers?
295  {
296  do
297  {
298  v1 = 2.0 * Next() - 1.0; // Pick two coordinates from
299  v2 = 2.0 * Next() - 1.0; // -1 to +1 and check if they
300  rsq = v1*v1 + v2*v2; // are in unit circle...
301  }
302  while (rsq >= 1.0 || rsq == 0.0);
303 
304  fac = sqrt(-2.0 * log(rsq)/rsq); // Apply the Box-Muller
305  normStore = v1 * fac; // transformation and save
306  normSaved = 1; // one deviate for next time
307  return v2 * fac;
308  }
309  else
310  {
311  normSaved = 0;
312  return normStore;
313  }
314 }
315 
316 void Random::Choose(int * array, int n, int k)
317 {
318  int choices = 1, others = 0;
319 
320  if (k > n / 2)
321  {
322  choices = 0;
323  others = 1;
324  k = n - k;
325  }
326 
327  for (int i = 0; i < n; i++)
328  array[i] = others;
329 
330  while (k > 0)
331  {
332  int i = NextInt() % n;
333 
334  if (array[i] == choices) continue;
335 
336  array[i] = choices;
337  k--;
338  }
339 }
340 
341 void Random::Choose(int * array, float * weights, int n, int k)
342 {
343  int choices = 1, others = 0;
344 
345  if (k > n / 2)
346  {
347  choices = 0;
348  others = 1;
349  k = n - k;
350  }
351 
352  // First calculate cumulative sums of weights ...
353  float * cumulative = new float [n + 1];
354 
355  cumulative[0] = 0;
356  for (int i = 1; i <= n; i++)
357  cumulative[i] = cumulative[i - 1] + weights[i - 1];
358 
359  float & sum = cumulative[n], reject = 0.0;
360 
361  for (int i = 0; i < n; i++)
362  array[i] = others;
363 
364  while (k > 0)
365  {
366  float weight = Next() * sum;
367 
368  int hi = n, lo = 0, i = 0;
369 
370  while (hi >= lo)
371  {
372  i = (hi + lo) / 2;
373 
374  if (cumulative[i + 1] <= weight)
375  lo = i + 1;
376  else if (cumulative[i] >= weight)
377  hi = i - 1;
378  else break;
379  }
380 
381  if (array[i] == choices) continue;
382 
383  array[i] = choices;
384  reject += weights[i];
385 
386  // After selecting a substantial number of elements, update the cumulative
387  // distribution -- to ensure that at least half of our samples produce a hit
388  if (reject > sum * 0.50)
389  {
390  cumulative[0] = 0;
391  for (int i = 1; i <= n; i++)
392  if (array[i] != choices)
393  cumulative[i] = cumulative[i - 1] + weights[i - 1];
394  else
395  cumulative[i] = cumulative[i - 1];
396 
397  reject = 0.0;
398  sum = cumulative[n];
399  }
400 
401  k--;
402  }
403 
404  delete [] cumulative;
405 }
406 
407 Random globalRandom;
408 
409 
Definition: Random.h:69