Commit b07bc1f0 authored by Andrew Heather's avatar Andrew Heather
Browse files

ENH: Random numbers - consolidated Random classes

- old Random class deprecated
- cachedRandom renamed Random
parent 88a03de2
......@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -25,173 +25,237 @@ License
#include "Random.H"
#include "OSspecific.H"
#include "PstreamReduceOps.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * //
#if INT_MAX != 2147483647
# error "INT_MAX != 2147483647"
# error "The random number generator may not work!"
#endif
Foam::scalar Foam::Random::scalar01()
{
return osRandomDouble(buffer_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::Random::Random(const label seed)
:
buffer_(osRandomBufferSize()),
seed_(seed),
hasGaussSample_(false),
gaussSample_(0)
{
if (seed > 1)
{
Seed = seed;
}
else
{
Seed = 1;
}
osRandomSeed(Seed);
// Initialise the random number generator
osRandomSeed(seed_, buffer_);
}
int Foam::Random::bit()
Foam::Random::Random(const Random& r, const bool reset)
:
buffer_(r.buffer_),
seed_(r.seed_),
hasGaussSample_(r.hasGaussSample_),
gaussSample_(r.gaussSample_)
{
if (osRandomInteger() > INT_MAX/2)
if (reset)
{
return 1;
}
else
{
return 0;
hasGaussSample_ = false;
gaussSample_ = 0;
// Re-initialise the samples
osRandomSeed(seed_, buffer_);
}
}
Foam::scalar Foam::Random::scalar01()
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::Random::~Random()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::Random::reset(const label seed)
{
return osRandomDouble();
seed_ = seed;
osRandomSeed(seed_, buffer_);
}
Foam::vector Foam::Random::vector01()
int Foam::Random::bit()
{
vector rndVec;
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
rndVec.component(cmpt) = scalar01();
}
return rndVec;
return osRandomInteger(buffer_) & 1;
}
Foam::sphericalTensor Foam::Random::sphericalTensor01()
template<>
Foam::scalar Foam::Random::sample01()
{
sphericalTensor rndTen;
rndTen.ii() = scalar01();
return scalar01();
}
return rndTen;
template<>
Foam::label Foam::Random::sample01()
{
return round(scalar01());
}
Foam::symmTensor Foam::Random::symmTensor01()
template<>
Foam::scalar Foam::Random::GaussNormal()
{
symmTensor rndTen;
for (direction cmpt=0; cmpt<symmTensor::nComponents; cmpt++)
if (hasGaussSample_)
{
rndTen.component(cmpt) = scalar01();
hasGaussSample_ = false;
return gaussSample_;
}
else
{
scalar rsq, v1, v2;
do
{
v1 = 2*scalar01() - 1;
v2 = 2*scalar01() - 1;
rsq = sqr(v1) + sqr(v2);
} while (rsq >= 1 || rsq == 0);
return rndTen;
scalar fac = sqrt(-2*log(rsq)/rsq);
gaussSample_ = v1*fac;
hasGaussSample_ = true;
return v2*fac;
}
}
Foam::tensor Foam::Random::tensor01()
template<>
Foam::label Foam::Random::GaussNormal()
{
tensor rndTen;
for (direction cmpt=0; cmpt<tensor::nComponents; cmpt++)
{
rndTen.component(cmpt) = scalar01();
}
return round(GaussNormal<scalar>());
}
return rndTen;
template<>
Foam::scalar Foam::Random::position
(
const scalar& start,
const scalar& end
)
{
return start + scalar01()*(end - start);
}
Foam::label Foam::Random::integer(const label lower, const label upper)
template<>
Foam::label Foam::Random::position(const label& start, const label& end)
{
return lower + (osRandomInteger() % (upper+1-lower));
return start + round(scalar01()*(end - start));
}
Foam::vector Foam::Random::position(const vector& start, const vector& end)
template<>
Foam::scalar Foam::Random::globalSample01()
{
vector rndVec(start);
scalar value = -GREAT;
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
if (Pstream::master())
{
rndVec.component(cmpt) +=
scalar01()*(end.component(cmpt) - start.component(cmpt));
value = scalar01();
}
return rndVec;
Pstream::scatter(value);
return value;
}
void Foam::Random::randomise(scalar& s)
template<>
Foam::label Foam::Random::globalSample01()
{
s = scalar01();
}
scalar value = -GREAT;
if (Pstream::master())
{
value = scalar01();
}
void Foam::Random::randomise(vector& v)
{
v = vector01();
Pstream::scatter(value);
return round(value);
}
void Foam::Random::randomise(sphericalTensor& st)
template<>
Foam::scalar Foam::Random::globalGaussNormal()
{
st = sphericalTensor01();
}
scalar value = -GREAT;
if (Pstream::master())
{
value = GaussNormal<scalar>();
}
void Foam::Random::randomise(symmTensor& st)
{
st = symmTensor01();
Pstream::scatter(value);
return value;
}
void Foam::Random::randomise(tensor& t)
template<>
Foam::label Foam::Random::globalGaussNormal()
{
t = tensor01();
scalar value = -GREAT;
if (Pstream::master())
{
value = GaussNormal<scalar>();
}
Pstream::scatter(value);
return round(value);
}
Foam::scalar Foam::Random::GaussNormal()
template<>
Foam::scalar Foam::Random::globalPosition
(
const scalar& start,
const scalar& end
)
{
static int iset = 0;
static scalar gset;
scalar fac, rsq, v1, v2;
scalar value = -GREAT;
if (iset == 0)
if (Pstream::master())
{
do
{
v1 = 2.0*scalar01() - 1.0;
v2 = 2.0*scalar01() - 1.0;
rsq = sqr(v1) + sqr(v2);
} while (rsq >= 1.0 || rsq == 0.0);
value = scalar01()*(end - start);
}
fac = sqrt(-2.0*log(rsq)/rsq);
gset = v1*fac;
iset = 1;
Pstream::scatter(value);
return v2*fac;
}
else
{
iset = 0;
return start + value;
}
return gset;
template<>
Foam::label Foam::Random::globalPosition
(
const label& start,
const label& end
)
{
label value = labelMin;
if (Pstream::master())
{
value = round(scalar01()*(end - start));
}
Pstream::scatter(value);
return start + value;
}
......
......@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -25,24 +25,29 @@ Class
Foam::Random
Description
Simple random number generator.
Random number generator.
SourceFiles
RandomI.H
Random.C
RandomTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef Random_H
#define Random_H
#include "vector.H"
#include "tensor.H"
#include "List.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class Random;
/*---------------------------------------------------------------------------*\
Class Random Declaration
\*---------------------------------------------------------------------------*/
......@@ -51,60 +56,165 @@ class Random
{
// Private data
label Seed;
//- Buffer used by re-entrant random number generator
List<char> buffer_;
//- Initial random number seed
label seed_;
public:
//- Indicator, which tells if there is a stored gaussian sample
bool hasGaussSample_;
//- Stored sample value
scalar gaussSample_;
// Private Member Functions
//- Returns the current sample
scalar scalar01();
public:
// Constructors
//- Construct given seed
Random(const label);
//- Construct given seed and sample count
Random(const label seed = 123456);
//- Copy constructor with optional reset of sampleI
Random(const Random& r, const bool reset = false);
// Destructor
~Random();
// Member functions
int bit();
// Access
//- Scalar [0..1] (so including 0,1)
scalar scalar01();
//- Return const access to the initial random number seed
inline label seed() const;
//- Reset the random number generator seed
void reset(const label seed);
// Evaluation
// Random numbers
//- Return a random bit
int bit();
//- Return a sample whose components lie in the range 0-1
template<class Type>
Type sample01();
//- Return a sample whose components are normally distributed
// with zero mean and unity variance N(0, 1)
template<class Type>
Type GaussNormal();
//- Return a sample between start and end
template<class Type>
Type position(const Type& start, const Type& end);
//- Vector with every component scalar01
vector vector01();
//- Randomise value in the range 0-1
template<class Type>
void randomise01(Type& value);
//- sphericalTensor with every component scalar01
sphericalTensor sphericalTensor01();
//- Shuffle the values in the list
template<class Type>
void shuffle(UList<Type>& values);
//- symmTensor with every component scalar01
symmTensor symmTensor01();
//- Tensor with every component scalar01
tensor tensor01();
// Global random numbers - consistent across all processors
//- Label [lower..upper]
label integer(const label lower, const label upper);
//- Return a sample whose components lie in the range 0-1
template<class Type>
Type globalSample01();
vector position(const vector&, const vector&);
//- Return a sample whose components are normally distributed
// with zero mean and unity variance N(0, 1)
template<class Type>
Type globalGaussNormal();
void randomise(scalar&);
void randomise(vector&);
void randomise(sphericalTensor&);
void randomise(symmTensor&);
void randomise(tensor&);
//- Return a sample between start and end
template<class Type>
Type globalPosition(const Type& start, const Type& end);
//- Return a normal Gaussian randon number
// with zero mean and unity variance N(0, 1)
scalar GaussNormal();
//- Randomise value in the range 0-1
template<class Type>
void globalRandomise01(Type& value);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Template specialisations
template<>
scalar Random::sample01<scalar>();
template<>
label Random::sample01<label>();
template<>
scalar Random::GaussNormal<scalar>();
template<>
label Random::GaussNormal<label>();
template<>
scalar Random::position<scalar>
(
const scalar& start,
const scalar& end
);
template<>
label Random::position<label>(const label& start, const label& end);
template<>
scalar Random::globalSample01<scalar>();
template<>
label Random::globalSample01<label>();
template<>
scalar Random::globalGaussNormal<scalar>();
template<>
label Random::globalGaussNormal<label>();
template<>
scalar Random::globalPosition<scalar>
(
const scalar& start,
const scalar& end
);
template<>
label Random::globalPosition<label>(const label& start, const label& end);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "RandomI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "RandomTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
......@@ -23,11 +23,11 @@ License
\*---------------------------------------------------------------------------*/
#include "cachedRandom.H"
#include "Random.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
inline Foam::label Foam::cachedRandom::seed() const
inline Foam::label Foam::Random::seed() const
{
return seed_;
}
......
......@@ -23,13 +23,13 @@ License
\*---------------------------------------------------------------------------*/
#include "cachedRandom.H"
#include "Random.H"
#include "Pstream.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Type Foam::cachedRandom::sample01()
Type Foam::Random::sample01()
{
Type value;
for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
......@@ -42,7 +42,7 @@ Type Foam::cachedRandom::sample01()
template<class Type>
Type Foam::cachedRandom::GaussNormal()
Type Foam::Random::GaussNormal()
{
Type value;
for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
......@@ -55,7 +55,7 @@ Type Foam::cachedRandom::GaussNormal()
template<class Type>