From b85d0b5cb7b5d75dc9e986cac91f9d45b12aa524 Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Mon, 9 Apr 2018 11:08:34 +0200 Subject: [PATCH] ENH: provide Rand48 as generator in the expected C++11 form - this removes an OS-specific dependency (eg, drand48_r is not POSIX) and allows easier use of other random number generators. The Rand48 generator has identical behaviour and period as the lrand48() library routine, but holds its own seed and state (which makes it re-entrant) and can be combined with other random distributions. However, when using the modified form to obtain scalar values they will not be identical to what drand48() yields. This is because drand48() uses the raw 48-bit values to directly set the mantissa of an IEEE double where as the newer distribution normalizes based on the 32-bit value. STYLE: simplify code in Random::shuffle and use Swap --- applications/test/Random/Make/files | 3 + applications/test/Random/Make/options | 3 + applications/test/Random/Test-Random.C | 183 ++++++++++++++++++ src/OSspecific/POSIX/POSIX.C | 86 -------- src/OpenFOAM/include/OSspecific.H | 24 --- .../primitives/random/Random/Rand48.H | 146 ++++++++++++++ .../primitives/random/Random/Random.C | 84 +++----- .../primitives/random/Random/Random.H | 129 ++++++------ .../primitives/random/Random/RandomI.H | 26 ++- .../random/Random/RandomTemplates.C | 21 +- 10 files changed, 455 insertions(+), 250 deletions(-) create mode 100644 applications/test/Random/Make/files create mode 100644 applications/test/Random/Make/options create mode 100644 applications/test/Random/Test-Random.C create mode 100644 src/OpenFOAM/primitives/random/Random/Rand48.H diff --git a/applications/test/Random/Make/files b/applications/test/Random/Make/files new file mode 100644 index 0000000000..ab0e0592c0 --- /dev/null +++ b/applications/test/Random/Make/files @@ -0,0 +1,3 @@ +Test-Random.C + +EXE = $(FOAM_USER_APPBIN)/Test-Random diff --git a/applications/test/Random/Make/options b/applications/test/Random/Make/options new file mode 100644 index 0000000000..4c3dd783cb --- /dev/null +++ b/applications/test/Random/Make/options @@ -0,0 +1,3 @@ +EXE_INC = + +EXE_LIBS = diff --git a/applications/test/Random/Test-Random.C b/applications/test/Random/Test-Random.C new file mode 100644 index 0000000000..b4bd2f2552 --- /dev/null +++ b/applications/test/Random/Test-Random.C @@ -0,0 +1,183 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2018 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Application + Test-Random + +Description + Simple test for sequence of random numbers + +\*---------------------------------------------------------------------------*/ + +#include "Rand48.H" +#include "Random.H" + +#include +#include +#include + +#define TEST_RAW_IEEE + +// Construct a positive double with the 48 random bits distributed over +// its fractional part so the resulting FP number is [0.0,1.0). +// +// As per glibc erand48() implementation + +#ifdef TEST_RAW_IEEE +#include +double randomFraction(const uint64_t bits) +{ + // 48-bit value + unsigned short int xsubi[3]; + xsubi[0] = (bits & 0xffff); + xsubi[1] = ((bits >> 16) & 0xffff); + xsubi[2] = ((bits >> 32) & 0xffff); + + union ieee754_double temp; + + temp.ieee.negative = 0; + temp.ieee.exponent = IEEE754_DOUBLE_BIAS; + temp.ieee.mantissa0 = (xsubi[2] << 4) | (xsubi[1] >> 12); + temp.ieee.mantissa1 = ((xsubi[1] & 0xfff) << 20) | (xsubi[0] << 4); + + // The lower 4 bits of mantissa1 are always 0. + return temp.d - 1.0; +} +#endif + + +using namespace Foam; + +// Output with cout instead of Info to keep values unsigned on output +using std::cout; +using std::setw; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + // Allow multiple passes to ensure reset is working properly + const label maxIter = 1; + + std::default_random_engine deflt(123456); + std::mt19937 mtwist(123456); + std::uniform_real_distribution uniform01; + + { + Rand48 rnd(123456); + + for (label iter=0; iter < maxIter; ++iter) + { + rnd.seed(123456); + ::srand48(123456); + mtwist.seed(123456); + + Info<< nl << "32-bit random with seed = 123456" << nl; + + cout<< setw(12) << "Rand48()" + << setw(12) << "lrand48()" + << setw(12) << "mtwister" + << setw(12) << "default" << nl; + + for (int i=0; i<25; i++) + { + cout<< setw(12) << rnd() + << setw(12) << long(::lrand48()) + << setw(12) << long(mtwist()) + << setw(12) << long(deflt()) << nl; + } + } + } + + { + Random rnd(123456); + Rand48 manual(123456); + mtwist.seed(123456); + deflt.seed(123456); + + // Two passes to ensure that reset is working properly + for (label iter=0; iter < maxIter; ++iter) + { + ::srand48(123456); + rnd.reset(123456); + manual.seed(123456); + mtwist.seed(123456); + deflt.seed(123456); + + cout<< nl << "Random (Rand48) with seed = " << rnd.seed() + << " interval [0,1000]" << nl; + + cout<< setw(12) << "Rand48()" + << setw(12) << "drand48()"; + + #ifdef TEST_RAW_IEEE + cout<< setw(12) << "manual"; + #endif + cout<< setw(12) << "mtwister"; + cout<< nl; + + for (int i=0; i<25; i++) + { + cout<< setw(12) << (rnd.sample01()*1000) + << setw(12) << (drand48()*1000); + + #ifdef TEST_RAW_IEEE + cout<< setw(12) << (randomFraction(manual.raw())*1000); + #endif + cout<< setw(12) << (uniform01(mtwist)*1000); + cout<< setw(12) << (uniform01(deflt)*1000); + cout<< nl; + } + } + } + + { + Rand48 rnd1(123456); + Rand48 rnd2(123456); + ::srand48(123456); + + rnd2.discard(10); + + cout<< nl << "Rand48 - test with offset of 10" << nl; + + cout<< setw(12) << "Rand48()" + << setw(12) << "offset-10" + << setw(12) << "lrand48()" + << nl; + + for (int i=0; i<25; i++) + { + cout<< setw(12) << (rnd1()) + << setw(12) << (rnd2()) + << setw(12) << long(::lrand48()) + << nl; + } + } + + cout<< nl << "Done." << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/src/OSspecific/POSIX/POSIX.C b/src/OSspecific/POSIX/POSIX.C index c37e911aee..b66699bc18 100644 --- a/src/OSspecific/POSIX/POSIX.C +++ b/src/OSspecific/POSIX/POSIX.C @@ -66,14 +66,6 @@ Description #include #endif -#ifdef USE_RANDOM - #include - #if INT_MAX != 2147483647 - #error "INT_MAX != 2147483647" - #error "The random number generator may not work!" - #endif -#endif - // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam @@ -1666,84 +1658,6 @@ Foam::fileNameList Foam::dlLoaded() return libs; } -Foam::label Foam::osRandomBufferSize() -{ - #ifdef USE_RANDOM - return sizeof(random_data); - #else - return sizeof(drand48_data); - #endif -} - - -void Foam::osRandomSeed(const label seed) -{ - #ifdef USE_RANDOM - srandom((unsigned int)seed); - #else - srand48(seed); - #endif -} - - -Foam::label Foam::osRandomInteger() -{ - #ifdef USE_RANDOM - return random(); - #else - return lrand48(); - #endif -} - - -Foam::scalar Foam::osRandomDouble() -{ - #ifdef USE_RANDOM - return (scalar)random()/INT_MAX; - #else - return drand48(); - #endif -} - - -void Foam::osRandomSeed(const label seed, List& buffer) -{ - #ifdef USE_RANDOM - srandom_r((unsigned int)seed, reinterpret_cast(buffer.begin())); - #else - srand48_r(seed, reinterpret_cast(buffer.begin())); - #endif -} - - -Foam::label Foam::osRandomInteger(List& buffer) -{ - #ifdef USE_RANDOM - int32_t result; - random_r(reinterpret_cast(buffer.begin()), &result); - return result; - #else - long result; - lrand48_r(reinterpret_cast(buffer.begin()), &result); - return result; - #endif -} - - -Foam::scalar Foam::osRandomDouble(List& buffer) -{ - #ifdef USE_RANDOM - int32_t result; - random_r(reinterpret_cast(buffer.begin()), &result); - return (scalar)result/INT_MAX; - #else - double result; - drand48_r(reinterpret_cast(buffer.begin()), &result); - return result; - #endif - -} - static Foam::DynamicList> threads_; static Foam::DynamicList> mutexes_; diff --git a/src/OpenFOAM/include/OSspecific.H b/src/OpenFOAM/include/OSspecific.H index bdd87752c4..c5077f01a2 100644 --- a/src/OpenFOAM/include/OSspecific.H +++ b/src/OpenFOAM/include/OSspecific.H @@ -256,30 +256,6 @@ bool dlSymFound(void* handle, const std::string& symbol); fileNameList dlLoaded(); -// Low level random numbers. Use Random class instead. - -//- Seed random number generator. -void osRandomSeed(const label seed); - -//- Return random integer (uniform distribution between 0 and 2^31) -label osRandomInteger(); - -//- Return random double precision (uniform distribution between 0 and 1) -scalar osRandomDouble(); - -//- Return the size of the buffer for re-entrant random number generator -label osRandomBufferSize(); - -//- Seed random number generator. -void osRandomSeed(const label seed, List& buffer); - -//- Return random integer (uniform distribution between 0 and 2^31) -label osRandomInteger(List& buffer); - -//- Return random double precision (uniform distribution between 0 and 1) -scalar osRandomDouble(List& buffer); - - // Thread handling //- Allocate a thread diff --git a/src/OpenFOAM/primitives/random/Random/Rand48.H b/src/OpenFOAM/primitives/random/Random/Rand48.H new file mode 100644 index 0000000000..fc0a2aec88 --- /dev/null +++ b/src/OpenFOAM/primitives/random/Random/Rand48.H @@ -0,0 +1,146 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2018 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::Rand48 + +Description + A pseudo random number generator using the linear congruential algorithm + with the following parameters: + \verbatim + Xn+1 = (aXn + c) mod m, where n >= 0 + + a = 0x5DEECE66D, + c = 0xB + m = 2**48 + \endverbatim + + It delivers results identical to the \c lrand48() function that is + available on some systems. + +\*---------------------------------------------------------------------------*/ + +#ifndef Rand48_H +#define Rand48_H + +#include +#include + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class Rand48 Declaration +\*---------------------------------------------------------------------------*/ + +class Rand48 +{ + //- It is a linear_congruential_engine + typedef std::linear_congruential_engine + < + uint64_t, + 0x5DEECE66D, + 0xB, + (uint64_t(1) << 48) + > engine; + + + // Private Data + + //- The calculation engine + engine engine_; + + + // Private Member Functions + + //- Convert integers + static constexpr uint64_t convert(uint32_t x) noexcept + { + return (static_cast(x) << 16) | 0x330e; + } + +public: + + //- The type of the generated random value. + typedef uint32_t result_type; + + //- The default seed + static constexpr result_type default_seed = 1u; + + + // Constructors + + //- Construct with specified or default seed + Rand48(uint32_t val = default_seed) + : + engine_(convert(val)) + {} + + + // Member Functions + + //- The smallest value that the generator can produce + static constexpr uint32_t min() { return 0; } + + //- The largest value that the generator can produce + static constexpr uint32_t max() { return 0x7FFFFFFF; } + + //- Reseed the generator + void seed(uint32_t val = default_seed) + { + engine_.seed(convert(val)); + } + + //- Advance the state of the generator by \c z notches. + void discard(unsigned long long z) + { + engine_.discard(z); + } + + //- Get the next random number in the sequence and return its raw + //- 48-bit value + uint64_t raw() + { + return engine_(); + } + + + // Member Operators + + //- Get the next random number in the sequence + uint32_t operator()() + { + return static_cast(engine_() >> 17); + } +}; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +#endif + +// ************************************************************************* // + diff --git a/src/OpenFOAM/primitives/random/Random/Random.C b/src/OpenFOAM/primitives/random/Random/Random.C index c17f5c2668..307fcb0eb1 100644 --- a/src/OpenFOAM/primitives/random/Random/Random.C +++ b/src/OpenFOAM/primitives/random/Random/Random.C @@ -2,8 +2,8 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2017-2018 OpenCFD Ltd. + \\/ M anipulation | ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -24,35 +24,25 @@ License \*---------------------------------------------------------------------------*/ #include "Random.H" -#include "OSspecific.H" #include "PstreamReduceOps.H" -// * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * // - -Foam::scalar Foam::Random::scalar01() -{ - return osRandomDouble(buffer_); -} - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::Random::Random(const label seed) +Foam::Random::Random(const label seedValue) : - buffer_(osRandomBufferSize()), - seed_(seed), + seed_(seedValue), + generator_(seed_), + uniform01_(), hasGaussSample_(false), gaussSample_(0) -{ - // Initialise the random number generator - osRandomSeed(seed_, buffer_); -} +{} Foam::Random::Random(const Random& r, const bool reset) : - buffer_(r.buffer_), seed_(r.seed_), + generator_(r.generator_), + uniform01_(), hasGaussSample_(r.hasGaussSample_), gaussSample_(r.gaussSample_) { @@ -60,34 +50,13 @@ Foam::Random::Random(const Random& r, const bool reset) { hasGaussSample_ = false; gaussSample_ = 0; - - // Re-initialise the samples - osRandomSeed(seed_, buffer_); + generator_.seed(seed_); } } -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::Random::~Random() -{} - - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void Foam::Random::reset(const label seed) -{ - seed_ = seed; - osRandomSeed(seed_, buffer_); -} - - -int Foam::Random::bit() -{ - return osRandomInteger(buffer_) & 1; -} - - template<> Foam::scalar Foam::Random::sample01() { @@ -110,23 +79,24 @@ Foam::scalar Foam::Random::GaussNormal() hasGaussSample_ = false; return gaussSample_; } - else + + // Gaussian random number as per Knuth/Marsaglia. + // Input: two uniform random numbers, output: two Gaussian random numbers. + // cache one of the values for the next call. + scalar rsq, v1, v2; + do { - scalar rsq, v1, v2; - do - { - v1 = 2*scalar01() - 1; - v2 = 2*scalar01() - 1; - rsq = sqr(v1) + sqr(v2); - } while (rsq >= 1 || rsq == 0); + v1 = 2*scalar01() - 1; + v2 = 2*scalar01() - 1; + rsq = sqr(v1) + sqr(v2); + } while (rsq >= 1 || rsq == 0); - scalar fac = sqrt(-2*log(rsq)/rsq); + const scalar fac = sqrt(-2*log(rsq)/rsq); - gaussSample_ = v1*fac; - hasGaussSample_ = true; + gaussSample_ = v1*fac; + hasGaussSample_ = true; - return v2*fac; - } + return v2*fac; } @@ -174,16 +144,16 @@ Foam::scalar Foam::Random::globalSample01() template<> Foam::label Foam::Random::globalSample01() { - scalar value = -GREAT; + label value = labelMin; if (Pstream::master()) { - value = scalar01(); + value = round(scalar01()); } Pstream::scatter(value); - return round(value); + return value; } diff --git a/src/OpenFOAM/primitives/random/Random/Random.H b/src/OpenFOAM/primitives/random/Random/Random.H index fb348fdafa..58d489b3e7 100644 --- a/src/OpenFOAM/primitives/random/Random/Random.H +++ b/src/OpenFOAM/primitives/random/Random/Random.H @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2017-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -27,7 +27,6 @@ Class Description Random number generator. - SourceFiles RandomI.H Random.C @@ -38,15 +37,18 @@ SourceFiles #ifndef Random_H #define Random_H -#include "List.H" +#include "Rand48.H" +#include "label.H" +#include "scalar.H" +#include // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { -// Forward declaration of classes -class Random; +// Forward declarations +template class UList; /*---------------------------------------------------------------------------*\ Class Random Declaration @@ -56,98 +58,99 @@ class Random { // Private data - //- Buffer used by re-entrant random number generator - List buffer_; - //- Initial random number seed label seed_; - //- Indicator, which tells if there is a stored gaussian sample + //- Random number generator on the int32 interval [0, 2^31) + Rand48 generator_; + + //- Uniform distribution on the scalar interval [0, 1] + std::uniform_real_distribution uniform01_; + + //- Is there a gaussian sample cached? bool hasGaussSample_; - //- Stored sample value + //- The cached gaussian sample value scalar gaussSample_; // Private Member Functions - //- Returns the current sample - scalar scalar01(); + //- A uniformly distributed floating-point random number [0,1] + inline scalar scalar01(); public: // Constructors - //- Construct given seed and sample count - Random(const label seed = 123456); + //- Construct with seed value + Random(const label seedValue = 123456); - //- Copy constructor with optional reset of sampleI + //- Copy construct with optional reset of seed Random(const Random& r, const bool reset = false); // Destructor - ~Random(); + ~Random() = default; - // Member functions + // Member Functions - // Access + // Access - //- Return const access to the initial random number seed - inline label seed() const; + //- The initial random number seed that was used + inline label seed() const; - //- Reset the random number generator seed - void reset(const label seed); + //- Reset the random number generator seed. + inline void reset(const label seedValue); - // Evaluation + // Random numbers - // Random numbers + //- Return a random bit + inline int bit(); - //- Return a random bit - int bit(); + //- Return a sample whose components lie in the range 0-1 + template + Type sample01(); - //- Return a sample whose components lie in the range 0-1 - template - Type sample01(); + //- Return a sample whose components are normally distributed + // with zero mean and unity variance N(0, 1) + template + Type GaussNormal(); - //- Return a sample whose components are normally distributed - // with zero mean and unity variance N(0, 1) - template - Type GaussNormal(); + //- Return a sample between start and end + template + Type position(const Type& start, const Type& end); - //- Return a sample between start and end - template - Type position(const Type& start, const Type& end); + //- Randomise value in the range 0-1 + template + void randomise01(Type& value); - //- Randomise value in the range 0-1 - template - void randomise01(Type& value); - - //- Shuffle the values in the list - template - void shuffle(UList& values); + //- Shuffle the values in the list + template + void shuffle(UList& values); - // Global random numbers - consistent across all processors + // Global random numbers - consistent across all processors - //- Return a sample whose components lie in the range 0-1 - template - Type globalSample01(); + //- Return a sample whose components lie in the range 0-1 + template + Type globalSample01(); - //- Return a sample whose components are normally distributed - // with zero mean and unity variance N(0, 1) - template - Type globalGaussNormal(); + //- Return a sample whose components are normally distributed + // with zero mean and unity variance N(0, 1) + template + Type globalGaussNormal(); - //- Return a sample between start and end - template - Type globalPosition(const Type& start, const Type& end); + //- Return a sample between start and end + template + Type globalPosition(const Type& start, const Type& end); - //- Randomise value in the range 0-1 - template - void globalRandomise01(Type& value); + //- Randomise value in the range 0-1 + template + void globalRandomise01(Type& value); }; @@ -167,11 +170,7 @@ template<> label Random::GaussNormal