Commit 9d25770b authored by Vedran Miletić's avatar Vedran Miletić

Replaced custom Bjarne-inspired RNG with PCG

PCG (pcg-random.org) is a family of simple fast space-efficient
statistically good algorithms for random number generation. Unlike many
general-purpose RNGs, they are also hard to predict.

The C++ implementation follows the conventions of other C++11-style
generators (WG21 N3551) and is available under the Apache License
version 2.0.
parent fea9832a
// Random number generator classes. Taken from "The C++ Programming Language",
// Third Edition, Bjarne Stroustrup, Addison-Wesley, 1997,p686.
#ifndef _RANDINT_H_
#define _RANDINT_H_
// uniform distribution in the interval [0,max]
class Randint {
unsigned long randx;
public:
Randint(long s = 0) { randx = s; }
void seed(long s) { randx = s; }
// DM 22 Feb 1999 - return current seed
long GetSeed() { return randx; }
// magic numbers chosen to use 31 bits of a 32-bit long:
int abs(int x) { return x & 0x7fffffff; }
static double fintmax() { return 2147483648.0; }
int draw() { return randx = randx * 1103515245 + 12345; }
double fdraw() { return abs(draw()) / fintmax(); }
int operator()() { return abs(draw()); }
};
// uniform distribution in the interval [0:n[
class Urand : public Randint {
int n;
public:
Urand(int nn) { n = nn; }
int operator()() {
int r = n * ((int)fdraw());
return (r == n) ? n - 1 : r;
}
};
#endif //_RANDINT_H_
......@@ -17,7 +17,9 @@
#ifndef _RBTRAND_H_
#define _RBTRAND_H_
#include "RandInt.h"
#include <pcg_random.hpp>
#include <random>
#include "RbtCoord.h"
class RbtRand {
......@@ -34,12 +36,12 @@ public:
// Seed the random number generator
RBTDLL_EXPORT void Seed(int seed = 0);
// Seed the random number generator from the system clock
void SeedFromClock();
// Seed the random number generator from the random device
void SeedFromRandomDevice();
// Returns current seed
RBTDLL_EXPORT int GetSeed();
// Get a random double between 0 and 1 (inlined)
double GetRandom01() { return m_rand.fdraw(); }
// Get a random double between 0 and 1
RBTDLL_EXPORT double GetRandom01();
// Get a random integer between 0 and nMax-1
int GetRandomInt(int nMax);
// Get a random unit vector distributed evenly over the surface of a sphere
......@@ -48,7 +50,7 @@ public:
double GetCauchyRandom(double, double);
private:
Randint m_rand; // Random number generator
pcg32 m_rng; // Random number generator
};
///////////////////////////////////////
......
......@@ -73,27 +73,29 @@ srcRbt = [
'import/simplex/src/NMSearch.cxx'
]
pcg_cpp_dep = dependency('pcg-cpp', fallback : ['pcg', 'pcg_cpp_dep'])
library_soversion = '0'
libRbt = shared_library(
'Rbt', srcRbt, soversion : library_soversion,
'Rbt', srcRbt, dependencies : pcg_cpp_dep, soversion : library_soversion,
version : library_soversion + '.' + meson.project_version(),
include_directories : incRbtAll
)
executable(
'rbdock', 'src/exe/rbdock.cxx', link_with : libRbt,
include_directories : incRbt
dependencies : pcg_cpp_dep, include_directories : incRbt
)
executable(
'rbcavity', 'src/exe/rbcavity.cxx', link_with : libRbt,
include_directories : incRbt
dependencies : pcg_cpp_dep, include_directories : incRbt
)
executable(
'rbrms', 'src/exe/rbrms.cxx', link_with : libRbt, include_directories : incRbt
)
executable(
'rbmoegrid', 'src/exe/rbmoegrid.cxx', link_with : libRbt,
include_directories : incRbt
dependencies : pcg_cpp_dep, include_directories : incRbt
)
executable(
'rblist', 'src/exe/rblist.cxx', link_with : libRbt,
......@@ -101,7 +103,7 @@ executable(
)
executable(
'rbcalcgrid', 'src/exe/rbcalcgrid.cxx', link_with : libRbt,
include_directories : incRbt
dependencies : pcg_cpp_dep, include_directories : incRbt
)
executable(
'rbconvgrid', 'src/exe/rbconvgrid.cxx', link_with : libRbt,
......@@ -138,8 +140,8 @@ if tests_opt
'build/test/RbtChromTest.cxx', 'build/test/SearchTest.cxx'
]
unit_test = executable(
'unit-test', srcTest, dependencies : gtest_dep, link_with : libRbt,
include_directories : [incTest, incRbt]
'unit-test', srcTest, dependencies : [pcg_cpp_dep, gtest_dep],
link_with : libRbt, include_directories : [incTest, incRbt]
)
test(
'unit-test', unit_test,
......
......@@ -25,7 +25,7 @@ int main(int argc, char *argv[]) {
if (argc > 1)
theRand.Seed(std::atoi(argv[2]));
else
theRand.SeedFromClock();
theRand.SeedFromRandomDevice();
std::cout << "Seed: " << theRand.GetSeed() << std::endl;
RbtGPGenome::SetStructure(56, 3, 1, 1, 7, 1, 200, 100);
ifstream desc("descnames", ios::in);
......
......@@ -25,7 +25,7 @@ int main(int argc, char *argv[]) {
if (argc > 1)
theRand.Seed(std::atoi(argv[1]));
else
theRand.SeedFromClock();
theRand.SeedFromRandomDevice();
std::cout << "Seed: " << theRand.GetSeed() << std::endl;
std::cout << "Output taken from: \n";
std::string strInputFile;
......
......@@ -77,7 +77,8 @@ void PrintUsage(void) {
std::cout << "\t\t-T <traceLevel> - controls output level for debugging (0 = "
"minimal, >0 = more verbose)"
<< std::endl;
std::cout << "\t\t-s <rndSeed> - random number seed (default=from sys clock)"
std::cout << "\t\t-s <rndSeed> - random number seed (default=from "
"std::random_device)"
<< std::endl;
}
......@@ -136,7 +137,7 @@ int main(int argc, char *argv[]) {
bool bNegIonise(false);
bool bImplH(true); // if true, read only polar hydrogens from SD file, else
// read all H's present
bool bSeed(false); // Random number seed (default = from system clock)
bool bSeed(false); // Random number seed (default = from std::random_device)
int nSeed(0);
bool bTrace(false);
int iTrace(0); // Trace level, for debugging
......@@ -508,8 +509,13 @@ int main(int argc, char *argv[]) {
if (spMdlFileSource->isDataFieldPresent("REG_Number"))
std::cout << "REG_Num:" << spMdlFileSource->GetDataValue("REG_Number")
<< std::endl;
std::cout << std::setw(30) << "RANDOM_NUMBER_SEED:" << theRand.GetSeed()
<< std::endl;
std::cout << std::setw(30) << "RNG seed:";
if (bSeed) {
std::cout << nSeed;
} else {
std::cout << "std::random_device";
}
std::cout << std::endl;
// Create and register the ligand model
RbtModelPtr spLigand = prmFactory.CreateLigand(spMdlFileSource);
......
......@@ -10,8 +10,6 @@
* http://rdock.sourceforge.net/
***********************************************************************/
#include <ctime> //Time functions for initialising the random number generator from the system clock
#include "RbtRand.h"
#include "Singleton.h"
......@@ -20,11 +18,11 @@
RbtRand::RbtRand() {
// Seed the random number generator
// Fixed seed in debug mode
// Seed from system clock in release mode
// Seed from the random device in release mode
#ifdef _DEBUG
Seed();
#else //_DEBUG
SeedFromClock();
SeedFromRandomDevice();
#endif //_DEBUG
_RBTOBJECTCOUNTER_CONSTR_("RbtRand");
}
......@@ -37,18 +35,27 @@ RbtRand::~RbtRand() { _RBTOBJECTCOUNTER_DESTR_("RbtRand"); }
// Public methods
// Seed the random number generator
void RbtRand::Seed(int seed) { m_rand.seed(seed); }
void RbtRand::Seed(int seed) { m_rng.seed(seed); }
// Seed the random number generator from the system clock
void RbtRand::SeedFromClock() { m_rand.seed(std::time(nullptr)); }
// Seed the random number generator from the random device
void RbtRand::SeedFromRandomDevice() {
pcg_extras::seed_seq_from<std::random_device> seed_source;
m_rng.seed(seed_source);
}
// Returns current seed
int RbtRand::GetSeed() { return m_rand.GetSeed(); }
// Get a random double between 0 and 1
double RbtRand::GetRandom01() {
std::uniform_real_distribution<> uniform_dist(0.0, 1.0);
return uniform_dist(m_rng);
}
// Get a random integer between 0 and nMax-1
int RbtRand::GetRandomInt(int nMax) {
int r = nMax * m_rand.fdraw();
return (r == nMax) ? nMax - 1 : r;
if (nMax == 0) {
return nMax;
}
std::uniform_int_distribution<> uniform_dist(0, nMax - 1);
return uniform_dist(m_rng);
}
// Get a random unit vector distributed evenly over the surface of a sphere
......@@ -63,28 +70,14 @@ RbtVector RbtRand::GetRandomUnitVector() {
// Get a random number from the Normal distribution (mean, variance)
double RbtRand::GetGaussianRandom(double mean, double variance) {
for (;;) {
double u1 = GetRandom01();
double u2 = GetRandom01();
double v1 = 2 * u1 - 1;
double v2 = 2 * u2 - 1;
double w = (v1 * v1) + (v2 * v2);
if (w <= 1) {
double y = std::sqrt((-2 * std::log(w)) / w);
double x1 = v1 * y;
return (x1 * std::sqrt(variance) + mean);
}
}
std::normal_distribution<> gaussian_dist{mean, variance};
return gaussian_dist(m_rng);
}
// Get a random number from the Cauchy distribution (mean, variance)
double RbtRand::GetCauchyRandom(double mean, double variance) {
double v1 = GetGaussianRandom(mean, variance);
double v2 = GetGaussianRandom(mean, variance);
if (std::fabs(v2) > 0.001)
return (v1 / v2);
else
return 0.0;
std::cauchy_distribution<> cauchy_dist{mean, variance};
return cauchy_dist(m_rng);
}
///////////////////////////////////////
......
[wrap-file]
directory = pcg-cpp-0.98.1
source_url = https://github.com/imneme/pcg-cpp/archive/v0.98.1.tar.gz
source_filename = pcg-cpp-v0.98.1.tar.gz
source_hash = 25bb22f076e8c346fa28c2267ae3564b12122f1f4ab2d2a08ad8909dcd6319fd
patch_url = https://wrapdb.mesonbuild.com/v1/projects/pcg/0.98.1/1/get_zip
patch_filename = pcg-0.98.1-1-wrap.zip
patch_hash = 99206a1e4e54b2f12b6177fa22d453dcdc612f67bc5d89760051c8bed509f783
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment