Skip to content
Permalink
Browse files

[FEATURE]: class for random transect sampling in analysis lib

  • Loading branch information
mhugent committed Jan 23, 2014
1 parent f4fc134 commit 5a193df48395f8102e7d08505ea699239b11b359
@@ -34,7 +34,9 @@ SET(QGIS_ANALYSIS_SRCS
raster/qgsrastercalcnode.cpp
raster/qgsrastercalculator.cpp
raster/qgsrastermatrix.cpp
vector/mersenne-twister.cpp
vector/qgsgeometryanalyzer.cpp
vector/qgstransectsample.cpp
vector/qgszonalstatistics.cpp
vector/qgsoverlayanalyzer.cpp

@@ -0,0 +1,226 @@
/*
* The Mersenne Twister pseudo-random number generator (PRNG)
*
* This is an implementation of fast PRNG called MT19937,
* meaning it has a period of 2^19937-1, which is a Mersenne
* prime.
*
* This PRNG is fast and suitable for non-cryptographic code.
* For instance, it would be perfect for Monte Carlo simulations,
* etc.
*
* For all the details on this algorithm, see the original paper:
* http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf
*
* Written by Christian Stigen Larsen
* http://csl.sublevel3.org
*
* Distributed under the modified BSD license.
*
* 2012-01-11
*/

#include <stdio.h>
#include <stdint.h>
#include <limits>
#include "mersenne-twister.h"

/*
* We have an array of 624 32-bit values, and there are
* 31 unused bits, so we have a seed value of
* 624*32-31 = 19937 bits.
*/
static const unsigned SIZE = 624;
static const unsigned PERIOD = 397;
static const unsigned DIFF = SIZE-PERIOD;

static uint32_t MT[SIZE];
static unsigned index = 0;

#define M32(x) (0x80000000 & x) // 32nd Most Significant Bit
#define L31(x) (0x7FFFFFFF & x) // 31 Least Significant Bits
#define ODD(x) (x & 1) // Check if number is odd

#define UINT32_MAX std::numeric_limits<uint32_t>::max()

static inline void generate_numbers()
{
/*
* Originally, we had one loop with i going from [0, SIZE) and
* two modulus operations:
*
* for ( register unsigned i=0; i<SIZE; ++i ) {
* register uint32_t y = M32(MT[i]) | L31(MT[(i+1) % SIZE]);
* MT[i] = MT[(i + PERIOD) % SIZE] ^ (y>>1);
* if ( ODD(y) ) MT[i] ^= 0x9908b0df;
* }
*
* For performance reasons, we've unrolled the loop three times, thus
* mitigating the need for any modulus operations.
*
* Anyway, it seems this trick is old hat:
* http://www.quadibloc.com/crypto/co4814.htm
*
*/

static const uint32_t MATRIX[2] = {0, 0x9908b0df};
register uint32_t y, i;

// i = [0 ... 226]
for ( i=0; i<DIFF; ++i ) {
/*
* We're doing 226 = 113*2, an even number of steps, so we can
* safely unroll one more step here for speed:
*/
y = M32(MT[i]) | L31(MT[i+1]);
MT[i] = MT[i+PERIOD] ^ (y>>1) ^ MATRIX[ODD(y)];

++i;
y = M32(MT[i]) | L31(MT[i+1]);
MT[i] = MT[i+PERIOD] ^ (y>>1) ^ MATRIX[ODD(y)];
}

#define UNROLL \
y = M32(MT[i]) | L31(MT[i+1]); \
MT[i] = MT[i-DIFF] ^ (y>>1) ^ MATRIX[ODD(y)]; \
++i;

// i = [227 ... 622]
for ( i=DIFF; i<(SIZE-1); ) {
/*
* 623-227 = 396 = 2*2*3*3*11, so we can unroll this loop in any
* number that evenly divides 396 (2, 4, 6, etc).
*/

// 11 times
UNROLL; UNROLL; UNROLL;
UNROLL; UNROLL; UNROLL;

UNROLL; UNROLL; UNROLL;
UNROLL; UNROLL;
}

// i = [623]
y = M32(MT[SIZE-1]) | L31(MT[SIZE-1]);
MT[SIZE-1] = MT[PERIOD-1] ^ (y>>1) ^ MATRIX[ODD(y)];
}

extern "C" void seed(uint32_t value)
{
/*
* The equation below is a linear congruential generator (LCG),
* one of the oldest known pseudo-random number generator
* algorithms, in the form X_(n+1) = = (a*X_n + c) (mod m).
*
* We've implicitly got m=32 (mask + word size of 32 bits), so
* there is no need to explicitly use modulus.
*
* What is interesting is the multiplier a. The one we have
* below is 0x6c07865 --- 1812433253 in decimal, and is called
* the Borosh-Niederreiter multiplier for modulus 2^32.
*
* It is mentioned in passing in Knuth's THE ART OF COMPUTER
* PROGRAMMING, Volume 2, page 106, Table 1, line 13. LCGs are
* treated in the same book, pp. 10-26
*
* You can read the original paper by Borosh and Niederreiter
* as well. It's called OPTIMAL MULTIPLIERS FOR PSEUDO-RANDOM
* NUMBER GENERATION BY THE LINEAR CONGRUENTIAL METHOD (1983) at
* http://www.springerlink.com/content/n7765ku70w8857l7/
*
* You can read about LCGs at:
* http://en.wikipedia.org/wiki/Linear_congruential_generator
*
* From that page, it says:
* "A common Mersenne twister implementation, interestingly
* enough, uses an LCG to generate seed data.",
*
* Since our we're using 32-bits data types for our MT array,
* we can skip the masking with 0xFFFFFFFF below.
*/

MT[0] = value;
index = 0;

for ( register unsigned i=1; i<SIZE; ++i )
MT[i] = 0x6c078965*(MT[i-1] ^ MT[i-1]>>30) + i;
}

extern "C" uint32_t rand_u32()
{
if ( !index )
generate_numbers();

register uint32_t y = MT[index];

// Tempering
y ^= y>>11;
y ^= y<< 7 & 0x9d2c5680;
y ^= y<<15 & 0xefc60000;
y ^= y>>18;

if ( ++index == SIZE )
index = 0;

return y;
}

extern "C" int mt_rand()
{
/*
* PORTABILITY WARNING:
*
* rand_u32() uses all 32-bits for the pseudo-random number,
* but rand() must return a number from 0 ... RAND_MAX.
*
* We'll just assume that rand() only uses 31 bits worth of
* data, and that we're on a two's complement system.
*
* So, to output an integer compatible with rand(), we have
* two options: Either mask off the highest (32nd) bit, or
* shift right by one bit. Masking with 0x7FFFFFFF will be
* compatible with 64-bit MT[], so we'll just use that here.
*
*/
return static_cast<int>(0x7FFFFFFF & rand_u32());
}

extern "C" void mt_srand(unsigned value)
{
seed(static_cast<uint32_t>(value));
}

extern "C" float randf_cc()
{
return static_cast<float>(rand_u32())/UINT32_MAX;
}

extern "C" float randf_co()
{
return static_cast<float>(rand_u32())/(UINT32_MAX+1.0f);
}

extern "C" float randf_oo()
{
return (static_cast<float>(rand_u32())+0.5f)/(UINT32_MAX+1.0f);
}

extern "C" double randd_cc()
{
return static_cast<double>(rand_u32())/UINT32_MAX;
}

extern "C" double randd_co()
{
return static_cast<double>(rand_u32())/(UINT32_MAX+1.0);
}

extern "C" double randd_oo()
{
return (static_cast<double>(rand_u32())+0.5)/(UINT32_MAX+1.0);
}

extern "C" uint64_t rand_u64()
{
return static_cast<uint64_t>(rand_u32())<<32 | rand_u32();
}
@@ -0,0 +1,103 @@
/*
* The Mersenne Twister pseudo-random number generator (PRNG)
*
* This is an implementation of fast PRNG called MT19937,
* meaning it has a period of 2^19937-1, which is a Mersenne
* prime.
*
* This PRNG is fast and suitable for non-cryptographic code.
* For instance, it would be perfect for Monte Carlo simulations,
* etc.
*
* This code has been designed as a drop-in replacement for libc rand and
* srand(). If you need to mix them, you should encapsulate this code in a
* namespace.
*
* Written by Christian Stigen Larsen
* 2012-01-11 -- http://csl.sublevel3.org
*
* Distributed under the modified BSD license.
*/

#ifndef MERSENNE_TWISTER_H
#define MERSENNE_TWISTER_H

#include <stdint.h>

#ifdef __cplusplus
extern "C" {
#endif

/*
* Maximum number you can get from rand().
*/
#define RAND_MAX INT32_MAX

/*
* Initialize the number generator with given seed.
* (LIBC REPLACEMENT FUNCTION)
*/
void mt_srand(unsigned seed_value);

/*
* Extract a pseudo-random integer in the range 0 ... RAND_MAX.
* (LIBC REPLACEMENT FUNCTION)
*/
int mt_rand();

/*
* Extract a pseudo-random unsigned 32-bit integer in the range 0 ... UINT32_MAX
*/
uint32_t rand_u32();

/*
* Combine two unsigned 32-bit pseudo-random numbers into one 64-bit
*/
uint64_t rand_u64();

/*
* Initialize Mersenne Twister with given seed value.
*/
void seed(uint32_t seed_value);

/*
* Return a random float in the CLOSED range [0, 1]
* Mnemonic: randf_co = random float 0=closed 1=closed
*/
float randf_cc();

/*
* Return a random float in the OPEN range [0, 1>
* Mnemonic: randf_co = random float 0=closed 1=open
*/
float randf_co();

/*
* Return a random float in the OPEN range <0, 1>
* Mnemonic: randf_oo = random float 0=open 1=open
*/
float randf_oo();

/*
* Return a random double in the CLOSED range [0, 1]
* Mnemonic: randd_co = random double 0=closed 1=closed
*/
double randd_cc();

/*
* Return a random double in the OPEN range [0, 1>
* Mnemonic: randd_co = random double 0=closed 1=open
*/
double randd_co();

/*
* Return a random double in the OPEN range <0, 1>
* Mnemonic: randd_oo = random double 0=open 1=open
*/
double randd_oo();

#ifdef __cplusplus
} // extern "C"
#endif

#endif // MERSENNE_TWISTER_H

0 comments on commit 5a193df

Please sign in to comment.
You can’t perform that action at this time.