libStatGen Software 1
Loading...
Searching...
No Matches
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
89Random::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
101Random::~Random()
102{
103#ifndef __NO_MERSENNE
104 delete [] mt;
105#else
106 delete [] shuffler;
107#endif
108}
109
110void 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
133void 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
148int Random::Binary()
149{
150 return Next() > 0.5 ? 1 : 0;
151}
152
153#ifndef __NO_MERSENNE
154
155double 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
202unsigned 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
249double 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
270unsigned 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
290double 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
316void 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
341void 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
407Random globalRandom;
408
409