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