diff options
Diffstat (limited to 'contrib/lua-torch/torch7/lib/TH/THRandom.c')
-rw-r--r-- | contrib/lua-torch/torch7/lib/TH/THRandom.c | 272 |
1 files changed, 272 insertions, 0 deletions
diff --git a/contrib/lua-torch/torch7/lib/TH/THRandom.c b/contrib/lua-torch/torch7/lib/TH/THRandom.c new file mode 100644 index 000000000..86d721e7b --- /dev/null +++ b/contrib/lua-torch/torch7/lib/TH/THRandom.c @@ -0,0 +1,272 @@ +#include "THGeneral.h" +#include "THRandom.h" + +/* Code for the Mersenne Twister random generator.... */ +#define n _MERSENNE_STATE_N +#define m _MERSENNE_STATE_M + +/* Creates (unseeded) new generator*/ +static THGenerator* THGenerator_newUnseeded(void) +{ + THGenerator *self = THAlloc(sizeof(THGenerator)); + memset(self, 0, sizeof(THGenerator)); + self->left = 1; + self->seeded = 0; + self->normal_is_valid = 0; + return self; +} + +/* Creates new generator and makes sure it is seeded*/ +THGenerator* THGenerator_new(void) +{ + THGenerator *self = THGenerator_newUnseeded(); + THRandom_seed(self); + return self; +} + +THGenerator* THGenerator_copy(THGenerator *self, THGenerator *from) +{ + memcpy(self, from, sizeof(THGenerator)); + return self; +} + +void THGenerator_free(THGenerator *self) +{ + THFree(self); +} + +int THGenerator_isValid(THGenerator *_generator) +{ + if ((_generator->seeded == 1) && + (_generator->left > 0 && _generator->left <= n) && (_generator->next <= n)) + return 1; + + return 0; +} + +#ifndef _WIN32 +#include <sys/types.h> +#include <sys/stat.h> +#include <fcntl.h> +#include <unistd.h> + +static unsigned long readURandomLong() +{ + int randDev = open("/dev/urandom", O_RDONLY); + unsigned long randValue; + if (randDev < 0) { + THError("Unable to open /dev/urandom"); + } + ssize_t readBytes = read(randDev, &randValue, sizeof(randValue)); + if (readBytes < sizeof(randValue)) { + THError("Unable to read from /dev/urandom"); + } + close(randDev); + return randValue; +} +#endif // _WIN32 + +unsigned long THRandom_seed(THGenerator *_generator) +{ +#ifdef _WIN32 + unsigned long s = (unsigned long)time(0); +#else + unsigned long s = readURandomLong(); +#endif + THRandom_manualSeed(_generator, s); + return s; +} + +/* The next 4 methods are taken from http:www.math.keio.ac.jpmatumotoemt.html + Here is the copyright: + Some minor modifications have been made to adapt to "my" C... */ + +/* + A C-program for MT19937, with initialization improved 2002/2/10. + Coded by Takuji Nishimura and Makoto Matsumoto. + This is a faster version by taking Shawn Cokus's optimization, + Matthe Bellew's simplification, Isaku Wada's double version. + + Before using, initialize the state by using init_genrand(seed) + or init_by_array(init_key, key_length). + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + + Any feedback is very welcome. + http://www.math.keio.ac.jp/matumoto/emt.html + email: matumoto@math.keio.ac.jp +*/ + +/* Macros for the Mersenne Twister random generator... */ +/* Period parameters */ +/* #define n 624 */ +/* #define m 397 */ +#define MATRIX_A 0x9908b0dfUL /* constant vector a */ +#define UMASK 0x80000000UL /* most significant w-r bits */ +#define LMASK 0x7fffffffUL /* least significant r bits */ +#define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) ) +#define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL)) +/*********************************************************** That's it. */ + +void THRandom_manualSeed(THGenerator *_generator, unsigned long the_seed_) +{ + int j; + + /* This ensures reseeding resets all of the state (i.e. state for Gaussian numbers) */ + THGenerator *blank = THGenerator_newUnseeded(); + THGenerator_copy(_generator, blank); + THGenerator_free(blank); + + _generator->the_initial_seed = the_seed_; + _generator->state[0] = _generator->the_initial_seed & 0xffffffffUL; + for(j = 1; j < n; j++) + { + _generator->state[j] = (1812433253UL * (_generator->state[j-1] ^ (_generator->state[j-1] >> 30)) + j); + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + /* In the previous versions, mSBs of the seed affect */ + /* only mSBs of the array state[]. */ + /* 2002/01/09 modified by makoto matsumoto */ + _generator->state[j] &= 0xffffffffUL; /* for >32 bit machines */ + } + _generator->left = 1; + _generator->seeded = 1; +} + +unsigned long THRandom_initialSeed(THGenerator *_generator) +{ + return _generator->the_initial_seed; +} + +void THRandom_nextState(THGenerator *_generator) +{ + unsigned long *p = _generator->state; + int j; + + _generator->left = n; + _generator->next = 0; + + for(j = n-m+1; --j; p++) + *p = p[m] ^ TWIST(p[0], p[1]); + + for(j = m; --j; p++) + *p = p[m-n] ^ TWIST(p[0], p[1]); + + *p = p[m-n] ^ TWIST(p[0], _generator->state[0]); +} + +unsigned long THRandom_random(THGenerator *_generator) +{ + unsigned long y; + + if (--(_generator->left) == 0) + THRandom_nextState(_generator); + y = *(_generator->state + (_generator->next)++); + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680UL; + y ^= (y << 15) & 0xefc60000UL; + y ^= (y >> 18); + + return y; +} + +/* generates a random number on [0,1)-double-interval */ +static double __uniform__(THGenerator *_generator) +{ + /* divided by 2^32 */ + return (double)THRandom_random(_generator) * (1.0/4294967296.0); +} + +/********************************************************* + + Thanks *a lot* Takuji Nishimura and Makoto Matsumoto! + + Now my own code... + +*********************************************************/ + +double THRandom_uniform(THGenerator *_generator, double a, double b) +{ + return(__uniform__(_generator) * (b - a) + a); +} + +double THRandom_normal(THGenerator *_generator, double mean, double stdv) +{ + THArgCheck(stdv > 0, 2, "standard deviation must be strictly positive"); + + /* This is known as the Box-Muller method */ + if(!_generator->normal_is_valid) + { + _generator->normal_x = __uniform__(_generator); + _generator->normal_y = __uniform__(_generator); + _generator->normal_rho = sqrt(-2. * log(1.0-_generator->normal_y)); + _generator->normal_is_valid = 1; + } + else + _generator->normal_is_valid = 0; + + if(_generator->normal_is_valid) + return _generator->normal_rho*cos(2.*M_PI*_generator->normal_x)*stdv+mean; + else + return _generator->normal_rho*sin(2.*M_PI*_generator->normal_x)*stdv+mean; +} + +double THRandom_exponential(THGenerator *_generator, double lambda) +{ + return(-1. / lambda * log(1-__uniform__(_generator))); +} + +double THRandom_cauchy(THGenerator *_generator, double median, double sigma) +{ + return(median + sigma * tan(M_PI*(__uniform__(_generator)-0.5))); +} + +/* Faut etre malade pour utiliser ca. + M'enfin. */ +double THRandom_logNormal(THGenerator *_generator, double mean, double stdv) +{ + THArgCheck(stdv > 0, 2, "standard deviation must be strictly positive"); + return(exp(THRandom_normal(_generator, mean, stdv))); +} + +int THRandom_geometric(THGenerator *_generator, double p) +{ + THArgCheck(p > 0 && p < 1, 1, "must be > 0 and < 1"); + return((int)(log(1-__uniform__(_generator)) / log(p)) + 1); +} + +int THRandom_bernoulli(THGenerator *_generator, double p) +{ + THArgCheck(p >= 0 && p <= 1, 1, "must be >= 0 and <= 1"); + return(__uniform__(_generator) <= p); +} |