1 /*
2  * Copyright 2011, Ben Langmead <langmea@cs.jhu.edu>
3  *
4  * This file is part of Bowtie 2.
5  *
6  * Bowtie 2 is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * Bowtie 2 is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with Bowtie 2.  If not, see <http://www.gnu.org/licenses/>.
18  */
19 
20 #ifndef RANDOM_GEN_H_
21 #define RANDOM_GEN_H_
22 
23 #include <stdint.h>
24 #include "assert_helpers.h"
25 
26 //#define MERSENNE_TWISTER
27 
28 #ifndef MERSENNE_TWISTER
29 
30 /**
31  * Simple pseudo-random linear congruential generator, a la Numerical
32  * Recipes.
33  */
34 class RandomSource {
35 public:
36 	static const uint32_t DEFUALT_A = 1664525;
37 	static const uint32_t DEFUALT_C = 1013904223;
38 
RandomSource()39 	RandomSource() :
40 		a(DEFUALT_A), c(DEFUALT_C), inited_(false) { }
RandomSource(uint32_t _last)41 	RandomSource(uint32_t _last) :
42 		a(DEFUALT_A), c(DEFUALT_C), last(_last), inited_(true) { }
RandomSource(uint32_t _a,uint32_t _c)43 	RandomSource(uint32_t _a, uint32_t _c) :
44 		a(_a), c(_c), inited_(false) { }
45 
46 	void init(uint32_t seed = 0) {
47 		last = seed;
48 		inited_ = true;
49 		lastOff = 30;
50 	}
51 
nextU32()52 	uint32_t nextU32() {
53 		assert(inited_);
54 		uint32_t ret;
55 		last = a * last + c;
56 		ret = last >> 16;
57 		last = a * last + c;
58 		ret ^= last;
59 		lastOff = 0;
60 		return ret;
61 	}
62 
nextU64()63     uint64_t nextU64() {
64 		assert(inited_);
65 		uint64_t first = nextU32();
66 		first = first << 32;
67 		uint64_t second = nextU32();
68 		return first | second;
69 	}
70 
71 
nextSizeT()72 	size_t nextSizeT() {
73 		if(sizeof(size_t) == 4) {
74 			return nextU32();
75 		} else {
76 			return nextU64();
77 		}
78 	}
79 
80 	/**
81 	 * Return a pseudo-random unsigned 32-bit integer sampled uniformly
82 	 * from [lo, hi].
83 	 */
nextU32Range(uint32_t lo,uint32_t hi)84 	uint32_t nextU32Range(uint32_t lo, uint32_t hi) {
85 		uint32_t ret = lo;
86 		if(hi > lo) {
87 			ret += (nextU32() % (hi-lo+1));
88 		}
89 		return ret;
90 	}
91 
92 	/**
93 	 * Get next 2-bit unsigned integer.
94 	 */
nextU2()95 	uint32_t nextU2() {
96 		assert(inited_);
97 		if(lastOff > 30) {
98 			nextU32();
99 		}
100 		uint32_t ret = (last >> lastOff) & 3;
101 		lastOff += 2;
102 		return ret;
103 	}
104 
105 	/**
106 	 * Get next boolean.
107 	 */
nextBool()108 	bool nextBool() {
109 		assert(inited_);
110 		if(lastOff > 31) {
111 			nextU32();
112 		}
113 		uint32_t ret = (last >> lastOff) & 1;
114 		lastOff++;
115 		return ret;
116 	}
117 
118 	/**
119 	 * Return an unsigned int chosen by picking randomly from among
120 	 * options weighted by probabilies supplied as the elements of the
121 	 * 'weights' array of length 'numWeights'.  The weights should add
122 	 * to 1.
123 	 */
nextFromProbs(const float * weights,size_t numWeights)124 	uint32_t nextFromProbs(
125 		const float* weights,
126 		size_t numWeights)
127 	{
128 		float f = nextFloat();
129 		float tot = 0.0f; // total weight seen so far
130 		for(uint32_t i = 0; i < numWeights; i++) {
131 			tot += weights[i];
132 			if(f < tot) return i;
133 		}
134 		return (uint32_t)(numWeights-1);
135 	}
136 
nextFloat()137 	float nextFloat() {
138 		assert(inited_);
139 		return (float)nextU32() / (float)0xffffffff;
140 	}
141 
142 	static uint32_t nextU32(uint32_t last,
143 	                        uint32_t a = DEFUALT_A,
144 	                        uint32_t c = DEFUALT_C)
145 	{
146 		return (a * last) + c;
147 	}
148 
currentA()149 	uint32_t currentA() const { return a; }
currentC()150 	uint32_t currentC() const { return c; }
currentLast()151 	uint32_t currentLast() const { return last; }
152 
153 private:
154 	uint32_t a;
155 	uint32_t c;
156 	uint32_t last;
157 	uint32_t lastOff;
158 	bool inited_;
159 };
160 
161 #else
162 
163 class RandomSource { // Mersenne Twister random number generator
164 
165 public:
166 
167 	// default constructor: uses default seed only if this is the first instance
RandomSource()168 	RandomSource() {
169 		reset();
170 	}
171 
172 	// constructor with 32 bit int as seed
RandomSource(uint32_t s)173 	RandomSource(uint32_t s) {
174 		init(s);
175 	}
176 
177 	// constructor with array of size 32 bit ints as seed
RandomSource(const uint32_t * array,int size)178 	RandomSource(const uint32_t* array, int size) {
179 		init(array, size);
180 	}
181 
reset()182 	void reset() {
183 		state_[0] = 0;
184 		p_ = 0;
185 		inited_ = false;
186 	}
187 
~RandomSource()188 	virtual ~RandomSource() { }
189 
190 	// the two seed functions
191 	void init(uint32_t); // seed with 32 bit integer
192 	void init(const uint32_t*, int size); // seed with array
193 
194 	/**
195 	 * Return next 1-bit unsigned integer.
196 	 */
nextBool()197 	bool nextBool() {
198 		return (nextU32() & 1) == 0;
199 	}
200 
201 	/**
202 	 * Get next unsigned 32-bit integer.
203 	 */
nextU32()204 	inline uint32_t nextU32() {
205 		assert(inited_);
206 		if(p_ == n) {
207 			gen_state(); // new state vector needed
208 		}
209 		// gen_state() is split off to be non-inline, because it is only called once
210 		// in every 624 calls and otherwise irand() would become too big to get inlined
211 		uint32_t x = state_[p_++];
212 		x ^= (x >> 11);
213 		x ^= (x << 7) & 0x9D2C5680UL;
214 		x ^= (x << 15) & 0xEFC60000UL;
215 		x ^= (x >> 18);
216 		return x;
217 	}
218 
219 	/**
220 	 * Return next float between 0 and 1.
221 	 */
nextFloat()222 	float nextFloat() {
223 		assert(inited_);
224 		return (float)nextU32() / (float)0xffffffff;
225 	}
226 
227 protected: // used by derived classes, otherwise not accessible; use the ()-operator
228 
229 	static const int n = 624, m = 397; // compile time constants
230 
231 	// the variables below are static (no duplicates can exist)
232 	uint32_t state_[n]; // state vector array
233 	int p_; // position in state array
234 
235 	bool inited_; // true if init function has been called
236 
237 	// private functions used to generate the pseudo random numbers
twiddle(uint32_t u,uint32_t v)238 	uint32_t twiddle(uint32_t u, uint32_t v) {
239 		return (((u & 0x80000000UL) | (v & 0x7FFFFFFFUL)) >> 1) ^ ((v & 1UL) ? 0x9908B0DFUL : 0x0UL);
240 	}
241 
242 	void gen_state(); // generate new state
243 
244 };
245 
246 #endif
247 
248 #endif /*RANDOM_GEN_H_*/
249