diff --git a/doc/obabel.1 b/doc/obabel.1 index 1411ccca1..f576bec0e 100644 --- a/doc/obabel.1 +++ b/doc/obabel.1 @@ -520,6 +520,13 @@ For further specific information and options, use .Fl H Ar format-type , e.g., .Fl Hcml +.Sh ENVIRONMENT +.Bl -tag -width flag +.It Ev OB_RANDOM_SEED +Seed for pseudo random number generator; random seed is used if unset or null. +This feature is effective if Open Babel is built with \fI-DOB_USE_OBRANDOMMT\fP cpreprocessor flag +(will be automatically enabled when version 4 is released). +.El .Sh EXAMPLES Standard conversion: .Dl "obabel ethanol.xyz \-Oethanol.pdb" diff --git a/doc/obdistgen.1 b/doc/obdistgen.1 index 786acf813..409b0f566 100644 --- a/doc/obdistgen.1 +++ b/doc/obdistgen.1 @@ -9,6 +9,13 @@ .Ar filename .Sh DESCRIPTION Generates SDFs with distance geometry. +.Sh ENVIRONMENT +.Bl -tag -width flag +.It Ev OB_RANDOM_SEED +Seed for pseudo random number generator; random seed is used if unset or null. +This feature is effective if Open Babel is built with \fI-DOB_USE_OBRANDOMMT\fP cpreprocessor flag +(will be automatically enabled when version 4 is released). +.El .Sh EXAMPLES .Dl "obdistgen molecule.smi > molecule.sdf" .Pp diff --git a/doc/obgen.1 b/doc/obgen.1 index 7901865a5..9ba2d9a4c 100644 --- a/doc/obgen.1 +++ b/doc/obgen.1 @@ -24,6 +24,13 @@ will give all options including the available forcefields. .It Fl ff Ar forcefield Select the forcefield .El +.Sh ENVIRONMENT +.Bl -tag -width flag +.It Ev OB_RANDOM_SEED +Seed for pseudo random number generator; random seed is used if unset or null. +This feature is effective if Open Babel is built with \fI-DOB_USE_OBRANDOMMT\fP cpreprocessor flag +(will be automatically enabled when version 4 is released). +.El .Sh EXAMPLES .Pp View the possible options, including available forcefields: diff --git a/doc/obrotamer.1 b/doc/obrotamer.1 index adc47e4ba..b8f1182a1 100644 --- a/doc/obrotamer.1 +++ b/doc/obrotamer.1 @@ -17,6 +17,13 @@ rotamers are not conformers -- that is, does not perform geometry optimization after generating the rotamer structure. The obminimize tool can do geometry optimization using molecular mechanics. +.Sh ENVIRONMENT +.Bl -tag -width flag +.It Ev OB_RANDOM_SEED +Seed for pseudo random number generator; random seed is used if unset or null. +This feature is effective if Open Babel is built with \fI-DOB_USE_OBRANDOMMT\fP cpreprocessor flag +(will be automatically enabled when version 4 is released). +.El .Sh EXAMPLES .Dl "obrotamer baseconformer.sdf > rotamer1.sdf" .Pp diff --git a/include/openbabel/conformersearch.h b/include/openbabel/conformersearch.h index 65bb7ca09..733a8e382 100644 --- a/include/openbabel/conformersearch.h +++ b/include/openbabel/conformersearch.h @@ -24,8 +24,14 @@ GNU General Public License for more details. #include #include #include +#include namespace OpenBabel { +#if !OB_USE_OBRANDOMMT + class OBRandom; +#else + class OBRandomMT; +#endif typedef std::vector RotorKey; typedef std::vector RotorKeys; @@ -443,7 +449,11 @@ namespace OpenBabel { std::vector > dynamic_niches; //!< The dynamically found niches std::vector niche_map; //!< Procide the sharing niche index, given the key inddex - void *d; // Opaque pointer - currently for storing OBRandom* which may be removed in future +#if !OB_USE_OBRANDOMMT + std::unique_ptr d; +#else + std::unique_ptr d; +#endif bool use_sharing; //!< Whether to use sharing or not. double alpha_share; //!< The alpha parameter in sharing function double sigma_share; //!< The sigma parameter in sharing function diff --git a/include/openbabel/math/matrix3x3.h b/include/openbabel/math/matrix3x3.h index f50bbc3b9..59e50aa08 100644 --- a/include/openbabel/math/matrix3x3.h +++ b/include/openbabel/math/matrix3x3.h @@ -36,8 +36,6 @@ GNU General Public License for more details. namespace OpenBabel { - class OBRandom; // class introduction in rand.h - // class introduction in matrix3x3.cpp class OBAPI matrix3x3 { diff --git a/src/config.h.cmake b/src/config.h.cmake index b2977e4cc..b438dab1d 100644 --- a/src/config.h.cmake +++ b/src/config.h.cmake @@ -200,3 +200,8 @@ #define TIME_WITH_SYS_TIME 0 #endif #endif + +// use OBRandomMT if true, otherwise deprecated OBRandom +#ifndef OB_USE_OBRANDOMMT +#define OB_USE_OBRANDOMMT (OB_VERSION >= OB_VERSION_CHECK(4, 0, 0)) +#endif diff --git a/src/conformersearch.cpp b/src/conformersearch.cpp index bd26ac011..c917b73f3 100644 --- a/src/conformersearch.cpp +++ b/src/conformersearch.cpp @@ -310,11 +310,12 @@ namespace OpenBabel { p_crossover = 0.7; niche_mating = 0.7; local_opt_rate = 3; - // For the moment 'd' is an opaque pointer to an instance of OBRandom*. - // In future, it could be a pointer to a structure storing all of the - // private variables. - d = (void*)new OBRandom(); - ((OBRandom*)d)->TimeSeed(); +#if !OB_USE_OBRANDOMMT + d.reset(new OBRandom()); + d->Reset(); +#else + d.reset(new OBRandomMT()); +#endif m_logstream = &std::cout; // Default logging send to standard output // m_logstream = NULL; m_printrotors = false; // By default, do not print rotors but perform the conformer search @@ -323,7 +324,6 @@ namespace OpenBabel { OBConformerSearch::~OBConformerSearch() { - delete (OBRandom*)d; } @@ -389,8 +389,12 @@ namespace OpenBabel { } // create initial population +#if !OB_USE_OBRANDOMMT OBRandom generator; - generator.TimeSeed(); + generator.Reset(); +#else + OBRandomMT generator{}; +#endif RotorKey rotorKey(m_rotorList.Size() + 1, 0); // indexed from 1 if (IsGood(rotorKey)) @@ -407,8 +411,13 @@ namespace OpenBabel { OBRotorIterator ri; OBRotor *rotor = m_rotorList.BeginRotor(ri); for (unsigned int i = 1; i < m_rotorList.Size() + 1; ++i, rotor = m_rotorList.NextRotor(ri)) { - if (generator.NextInt() % m_mutability == 0) - rotorKey[i] = generator.NextInt() % rotor->GetResolution().size(); + if (generator.UniformInt(0, m_mutability - 1) == 0) { +#if !OB_USE_OBRANDOMMT + rotorKey[i] = generator.UniformInt(0, rotor->GetResolution().size() - 1u); +#else + rotorKey[i] = generator.UniformInt(0, rotor->GetResolution().size() - 1u); +#endif + } } // duplicates are always rejected if (!IsUniqueKey(m_rotorKeys, rotorKey)) @@ -455,8 +464,12 @@ namespace OpenBabel { void OBConformerSearch::NextGeneration() { // create next generation population +#if !OB_USE_OBRANDOMMT OBRandom generator; - generator.TimeSeed(); + generator.Reset(); +#else + OBRandomMT generator{}; +#endif // generate the children int numConformers = m_rotorKeys.size(); @@ -473,8 +486,13 @@ namespace OpenBabel { OBRotorIterator ri; OBRotor *rotor = m_rotorList.BeginRotor(ri); for (unsigned int i = 1; i < m_rotorList.Size() + 1; ++i, rotor = m_rotorList.NextRotor(ri)) { - if (generator.NextInt() % m_mutability == 0) - rotorKey[i] = generator.NextInt() % rotor->GetResolution().size(); // permutate gene + if (generator.UniformInt(0, m_mutability - 1) == 0) { +#if !OB_USE_OBRANDOMMT + rotorKey[i] = generator.UniformInt(0, rotor->GetResolution().size() - 1u); // permutate gene +#else + rotorKey[i] = generator.UniformInt(0, rotor->GetResolution().size() - 1u); // permutate gene +#endif + } } // duplicates are always rejected if (!IsUniqueKey(m_rotorKeys, rotorKey)) @@ -792,9 +810,13 @@ namespace OpenBabel { for (i = 1; i <= m_rotorList.Size(); ++i, rotor = m_rotorList.NextRotor(ri)) { neighbor = best; - new_val = ((OBRandom*)d)->NextInt() % rotor->GetResolution().size(); - while (new_val == best[i]) - new_val = ((OBRandom*)d)->NextInt() % rotor->GetResolution().size(); + do { +#if !OB_USE_OBRANDOMMT + new_val = d->UniformInt(0, rotor->GetResolution().size() - 1u); +#else + new_val = d->UniformInt(0, rotor->GetResolution().size() - 1u); +#endif + } while (new_val == best[i]); neighbor[i] = new_val; if (IsUniqueKey(backup_population, neighbor) && IsGood(neighbor)) m_rotorKeys.push_back (neighbor); @@ -846,30 +868,30 @@ namespace OpenBabel { return 0; // Make a 2-tournament selection to choose first parent - i = ((OBRandom*)d)->NextInt() % pop_size; - j = ((OBRandom*)d)->NextInt() % pop_size; + i = d->UniformInt(0u, pop_size - 1u); + j = d->UniformInt(0u, pop_size - 1u); parent1 = vshared_fitnes[i] > vshared_fitnes[j] ? i : j; iniche = niche_map[parent1]; if (iniche > -1) nsize = dynamic_niches[iniche].size (); // Belongs to a specific niche // Do we apply crossover here? - flag_crossover = (((OBRandom*)d)->NextFloat () <= p_crossover); - if (flag_crossover && (((OBRandom*)d)->NextFloat () <= niche_mating) && nsize > 1) + flag_crossover = d->Bernoulli(p_crossover); + if (flag_crossover && d->Bernoulli(niche_mating) && nsize > 1) { // Apply niche mating: draw second parent in the same niche, if its has // at least 2 members. Make a 2-tournament selection whithin this niche - rnd1 = ((OBRandom*)d)->NextInt() % nsize; + rnd1 = d->UniformInt(0u, nsize - 1u); i = dynamic_niches[iniche][rnd1]; - rnd2 = ((OBRandom*)d)->NextInt() % nsize; + rnd2 = d->UniformInt(0u, nsize - 1u); j = dynamic_niches[iniche][rnd2]; parent2 = vshared_fitnes[i] > vshared_fitnes[j] ? i : j; } else { // Draw second in the whole population - i = ((OBRandom*)d)->NextInt() % pop_size; - j = ((OBRandom*)d)->NextInt() % pop_size; + i = d->UniformInt(0u, pop_size - 1u); + j = d->UniformInt(0u, pop_size - 1u); parent2 = vshared_fitnes[i] > vshared_fitnes[j] ? i : j; } @@ -878,7 +900,7 @@ namespace OpenBabel { // Cross the 2 vectors: toss a coin for each position (i.e. uniform crossover) for (i = 1; i < key1.size(); i++) { - if (((OBRandom*)d)->NextInt() % 2) + if (d->UniformInt(0, 1) == 0) { // Copy parent1 to offspring 1 key1[i] = m_rotorKeys[parent1][i]; key2[i] = m_rotorKeys[parent2][i]; @@ -900,10 +922,20 @@ namespace OpenBabel { rotor = m_rotorList.BeginRotor(ri); for (i = 1; i <= m_rotorList.Size(); ++i, rotor = m_rotorList.NextRotor(ri)) { - if (((OBRandom*)d)->NextInt() % m_mutability == 0) - key1[i] = ((OBRandom*)d)->NextInt() % rotor->GetResolution().size(); - if (((OBRandom*)d)->NextInt() % m_mutability == 0) - key2[i] = ((OBRandom*)d)->NextInt() % rotor->GetResolution().size(); + if (d->UniformInt(0, m_mutability - 1) == 0) { +#if !OB_USE_OBRANDOMMT + key1[i] = d->UniformInt(0, rotor->GetResolution().size() - 1u); +#else + key1[i] = d->UniformInt(0, rotor->GetResolution().size() - 1u); +#endif + } + if (d->UniformInt(0, m_mutability - 1) == 0) { +#if !OB_USE_OBRANDOMMT + key2[i] = d->UniformInt(0, rotor->GetResolution().size() - 1u); +#else + key2[i] = d->UniformInt(0, rotor->GetResolution().size() - 1u); +#endif + } } if (IsUniqueKey(m_rotorKeys, key1) && IsGood(key1)) ret_code += 1; diff --git a/src/distgeom.cpp b/src/distgeom.cpp index 81fce1c19..4c2bd0191 100644 --- a/src/distgeom.cpp +++ b/src/distgeom.cpp @@ -1015,13 +1015,17 @@ namespace OpenBabel { unsigned int N = _mol.NumAtoms(); // random distance matrix Eigen::MatrixXd distMat = Eigen::MatrixXd::Zero(N, N); +#if !OB_USE_OBRANDOMMT OBRandom generator; - generator.TimeSeed(); + generator.Reset(); +#else + OBRandomMT generator{}; +#endif for (size_t i=0; iGetLowerBounds(i, j); double ub = _d->GetUpperBounds(i, j); - double v = generator.NextFloat() * (ub - lb) + lb; + double v = generator.UniformReal(lb, ub); distMat(i, j) = v; distMat(j, i) = v; } @@ -1174,8 +1178,12 @@ namespace OpenBabel { _mol.AddConformer(confCoord); _mol.SetConformer(_mol.NumConformers()); +#if !OB_USE_OBRANDOMMT OBRandom generator(true); // Use system rand() functions - generator.TimeSeed(); + generator.Reset(); +#else + OBRandomMT generator{}; +#endif if (_d->debug) { cerr << " max box size: " << _d->maxBoxSize << endl; diff --git a/src/forcefield.cpp b/src/forcefield.cpp index e314f002a..889fd29fb 100644 --- a/src/forcefield.cpp +++ b/src/forcefield.cpp @@ -1433,8 +1433,12 @@ namespace OpenBabel OBRotorIterator ri; OBRotor *rotor; +#if !OB_USE_OBRANDOMMT OBRandom generator; - generator.TimeSeed(); + generator.Reset(); +#else + OBRandomMT generator{}; +#endif _origLogLevel = _loglvl; if (_mol.GetCoordinates() == nullptr) @@ -1479,7 +1483,11 @@ namespace OpenBabel rotor = rl.BeginRotor(ri); for (unsigned int i = 1; i < rl.Size() + 1; ++i, rotor = rl.NextRotor(ri)) { // foreach rotor - rotorKey[i] = generator.NextInt() % rotor->GetResolution().size(); +#if !OB_USE_OBRANDOMMT + rotorKey[i] = generator.UniformInt(0, rotor->GetResolution().size() - 1u); +#else + rotorKey[i] = generator.UniformInt(0, rotor->GetResolution().size() - 1u); +#endif } rotamers.AddRotamer(rotorKey); } @@ -1603,8 +1611,12 @@ namespace OpenBabel OBRotorIterator ri; OBRotor *rotor; +#if !OB_USE_OBRANDOMMT OBRandom generator; - generator.TimeSeed(); + generator.Reset(); +#else + OBRandomMT generator{}; +#endif int origLogLevel = _loglvl; if (_mol.GetCoordinates() == nullptr) @@ -1752,11 +1764,10 @@ namespace OpenBabel for (unsigned int i = 1; i < rl.Size() + 1; ++i, rotor = rl.NextRotor(ri)) { // foreach rotor rotorKey[i] = -1; // default = don't change dihedral - randFloat = generator.NextFloat(); - if (randFloat < defaultRotor) // should we just leave this rotor with default setting? + if (generator.Bernoulli(defaultRotor)) // should we just leave this rotor with default setting? continue; - randFloat = generator.NextFloat(); + randFloat = generator.UniformReal(0.0, 1.0); total = 0.0; for (unsigned int j = 0; j < rotor->GetResolution().size(); j++) { if (randFloat > total && randFloat < (total+ rotorWeights[i][j])) { @@ -3416,11 +3427,14 @@ namespace OpenBabel void OBForceField::GenerateVelocities() { cout << "OBForceField::GenerateVelocities()" << endl; +#if !OB_USE_OBRANDOMMT OBRandom generator; - generator.TimeSeed(); + generator.Reset(); +#else + OBRandomMT generator{}; +#endif _ncoords = _mol.NumAtoms() * 3; int velocityIdx; - double velocity; _velocityPtr = new double[_ncoords]; memset(_velocityPtr, '\0', sizeof(double)*_ncoords); @@ -3429,32 +3443,20 @@ namespace OpenBabel if (!_constraints.IsFixed(a->GetIdx()) || (_fixAtom == a->GetIdx()) || (_ignoreAtom == a->GetIdx())) { velocityIdx = (a->GetIdx() - 1) * 3; - // add twelve random numbers between 0.0 and 1.0, - // subtract 6.0 from their sum, multiply with sqrt(kT/m) + // sqrt(kT/m) + const double sigma = sqrt((GAS_CONSTANT * _temp)/ (1000 * a->GetAtomicMass())); if (!_constraints.IsXFixed(a->GetIdx())) { - velocity = 0.0; - for (int i=0; i < 12; ++i) - velocity += generator.NextFloat(); - velocity -= 6.0; - velocity *= sqrt((GAS_CONSTANT * _temp)/ (1000 * a->GetAtomicMass())); + double velocity = generator.Normal(0.0, sigma); _velocityPtr[velocityIdx] = velocity; // x10: gromacs uses nm instead of A } if (!_constraints.IsYFixed(a->GetIdx())) { - velocity = 0.0; - for (int i=0; i < 12; ++i) - velocity += generator.NextFloat(); - velocity -= 6.0; - velocity *= sqrt((GAS_CONSTANT * _temp)/ (1000 * a->GetAtomicMass())); + double velocity = generator.Normal(0.0, sigma); _velocityPtr[velocityIdx+1] = velocity; // idem } if (!_constraints.IsZFixed(a->GetIdx())) { - velocity = 0.0; - for (int i=0; i < 12; ++i) - velocity += generator.NextFloat(); - velocity -= 6.0; - velocity *= sqrt((GAS_CONSTANT * _temp)/ (1000 * a->GetAtomicMass())); + double velocity = generator.Normal(0.0, sigma); _velocityPtr[velocityIdx+2] = velocity; // idem } } diff --git a/src/math/vector3.cpp b/src/math/vector3.cpp index a0f79c762..3fe1343b2 100644 --- a/src/math/vector3.cpp +++ b/src/math/vector3.cpp @@ -27,6 +27,10 @@ GNU General Public License for more details. #include "../rand.h" #include +#ifndef OB_USE_IMPROVED_RANDOM_UNIT_VECTOR +#define OB_USE_IMPROVED_RANDOM_UNIT_VECTOR (OB_VERSION >= OB_VERSION_CHECK(4, 0, 0)) +#endif + using namespace std; namespace OpenBabel @@ -84,26 +88,33 @@ namespace OpenBabel } /*! Replaces *this with a random unit vector, which is (supposed - to be) uniformly distributed over the unit sphere. Uses the - system number generator with a time seed. - + to be) uniformly distributed over the unit sphere. */ void vector3::randomUnitVector() { - OBRandom *ptr; +#if !OB_USE_OBRANDOMMT static OBRandom singleRand(true); - ptr = &singleRand; +#else + static OBRandomMT singleRand{}; +#endif +#if !OB_USE_IMPROVED_RANDOM_UNIT_VECTOR // obtain a random vector with 0.001 <= length^2 <= 1.0, normalize // the vector to obtain a random vector of length 1.0. double l; do { - this->Set(ptr->NextFloat()-0.5, ptr->NextFloat()-0.5, ptr->NextFloat()-0.5); + this->Set(singleRand.UniformReal(-0.5, 0.5), singleRand.UniformReal(-0.5, 0.5), singleRand.UniformReal(-0.5, 0.5)); l = length_2(); } while ( (l > 1.0) || (l < 1e-4) ); this->normalize(); +#else + double z = singleRand.UniformReal(-1.0, 1.0); + double phi = singleRand.UniformReal(0.0, 2 * M_PI); + double rho = sqrt(1.0 - z * z); + this->Set(rho * cos(phi), rho * sin(phi), z); +#endif } OBAPI ostream& operator<< ( ostream& co, const vector3& v ) diff --git a/src/rand.cpp b/src/rand.cpp index 0975cbf0f..5ce264c9a 100644 --- a/src/rand.cpp +++ b/src/rand.cpp @@ -27,6 +27,7 @@ GNU General Public License for more details. #include #include #include +#include #include "rand.h" @@ -520,6 +521,52 @@ namespace OpenBabel #endif } + double OBRandom::Normal(double mu, double sigma) { + const int n = 12; + double z = -0.5 * n; + for (int _ = 0; _ < n; _++) { + z += NextFloat(); + } + return mu + z * sigma; + } + + bool OBRandom::Bernoulli(double p) { + return UniformReal(0.0, 1.0) <= p; + } + + uint_fast64_t OBRandomMT::randomSeed() const { + auto ob_random_seed = std::getenv("OB_RANDOM_SEED"); + if (ob_random_seed && ob_random_seed[0] != '\0') { + return std::atol(ob_random_seed); + } +#if defined(WIN32) || defined(__MINGW32__) + // for VC++ do it this way + time_t ltime; + time(<ime); + return ltime; +#else + return std::random_device{}(); +#endif + } + + int OBRandomMT::NextInt() const { + return UniformInt(0, std::numeric_limits::max()); + } + + double OBRandomMT::NextFloat() const { + return UniformReal(0.0, 1.0); + } + + double OBRandomMT::Normal(double mu, double sigma) const { + std::normal_distribution n{mu, sigma}; + return n(*prng); + } + + bool OBRandomMT::Bernoulli(double p) const { + std::bernoulli_distribution b{p}; + return b(*prng); + } + } //end namespace OpenBabel //! \file rand.cpp diff --git a/src/rand.h b/src/rand.h index 2cb508ed6..9c0f83c3b 100644 --- a/src/rand.h +++ b/src/rand.h @@ -20,6 +20,8 @@ GNU General Public License for more details. #ifndef RAND_H #define RAND_H +#include +#include #include namespace OpenBabel @@ -62,7 +64,7 @@ namespace OpenBabel generator.Seed(10246);// Use a specific initial seed value for reproducing sequence \endcode **/ - class OBRandom + class OB_DEPRECATED_MSG("use OBRandomMT instead") OBRandom { DoubleType d; unsigned int m,a,c; @@ -83,10 +85,78 @@ namespace OpenBabel //! If sranddev is available (e.g., Mac OS X, BSD...) use this instead //! for more random seeds void TimeSeed(); + void Reset() { TimeSeed(); } //! \return a random integer int NextInt(); //! \return a random floating-point number between 0.0 and 1.0 double NextFloat(); + //! \param a lower bound + //! \param b upper bound + //! \warning this implementation has mudulo bias + int UniformInt(int a, int b) { return a + NextInt() % (b - a + 1); } + //! \param a lower bound + //! \param b upper bound + double UniformReal(double a, double b) { return a + (b - a) * NextFloat(); } + //! \param mu mean + //! \param sigma standard deviation + //! \warning normal distribution is approximated by Irwin–Hall distribution + double Normal(double mu, double sigma); + //! \param p probability of true + bool Bernoulli(double p = 0.5); + }; + + //! @class PRNG based on Mersenne Twister + class OBRandomMT { + public: + //! initialize by a random seed + OBRandomMT() : OBRandomMT{randomSeed()} {} + OB_DEPRECATED explicit OBRandomMT(bool) : OBRandomMT{} {} + //! @param seed initializes pseudo random number generator + explicit OBRandomMT(uint_fast64_t seed) { + prng.reset(new std::mt19937_64{seed}); + } + ~OBRandomMT() = default; + + //! @param seed reinitializes pseudo random number generator + OB_DEPRECATED_MSG("use Reset instead") void Seed(uint_fast64_t seed) { Reset(seed); } + OB_DEPRECATED_MSG("use Reset instead") void TimeSeed() { Reset(); } + //! reset underlying pseudo random number generator by random seed + void Reset() { Reset(randomSeed()); } + //! @param seed reinitializes pseudo random number generator + void Reset(uint_fast64_t seed) { prng.reset(new std::mt19937_64{seed}); } + //! @return a random integer in [0, 2^31 - 1] + OB_DEPRECATED_MSG("use UniformInt instead") int NextInt() const; + //! @return a random floating-point number in [0.0, 1.0) + OB_DEPRECATED_MSG("use UniformReal instead") double NextFloat() const; + //! @param a lower bound + //! @param b upper bound + //! @return a random integer in [a, b] + template + T UniformInt(T a, T b) const { + std::uniform_int_distribution u{a, b}; + return u(*prng); + } + //! @param a lower bound + //! @param b upper bound + //! @return a random floating-point number in [a, b) + template + T UniformReal(T a, T b) const { + std::uniform_real_distribution u{a, b}; + return u(*prng); + } + //! @param mu mean + //! @param sigma standard deviation + double Normal(double mu, double sigma) const; + //! @param p probability of true + bool Bernoulli(double p = 0.5) const; + + private: + //! @return (i) OB_RANDOM_SEED environment variable if set and non-null, + //! (ii) random number obtained by std::random_device if available, + //! (iii) otherwise the current time + uint_fast64_t randomSeed() const; + + std::unique_ptr prng; }; } // end namespace OpenBabel diff --git a/test/math.cpp b/test/math.cpp index 4d776431f..bbe815303 100644 --- a/test/math.cpp +++ b/test/math.cpp @@ -51,7 +51,11 @@ GNU General Public License for more details. using namespace std; using namespace OpenBabel; +#if !OB_USE_OBRANDOMMT OBRandom randomizer; +#else +OBRandomMT randomizer{}; +#endif int testCount = 0; int failedCount = 0; char currentFunc [256]; @@ -72,7 +76,7 @@ void verify_not_ok( const char *expr, int line ) void pickRandom( double & d ) { - d = randomizer.NextFloat() * 2.0 - 1.0; + d = randomizer.UniformReal(-1.0, 1.0); } void pickRandom( vector3 & v ) @@ -323,9 +327,9 @@ void testEigenvalues() for(int i=0; i<3; i++) for(int j=0; j<3; j++) Diagonal.Set(i, j, 0.0); - Diagonal.Set(0, 0, randomizer.NextFloat()); - Diagonal.Set(1, 1, Diagonal.Get(0,0)+fabs(randomizer.NextFloat())); - Diagonal.Set(2, 2, Diagonal.Get(1,1)+fabs(randomizer.NextFloat())); + Diagonal.Set(0, 0, randomizer.UniformReal(0.0, 1.0)); + Diagonal.Set(1, 1, Diagonal.Get(0,0) + randomizer.UniformReal(0.0, 1.0)); + Diagonal.Set(2, 2, Diagonal.Get(1,1) + randomizer.UniformReal(0.0, 1.0)); // test the isDiagonal() method VERIFY( Diagonal.isDiagonal() ); @@ -397,7 +401,7 @@ int math(int argc, char* argv[]) cout << "# math: repeating each test " << REPEAT << " times" << endl; - randomizer.TimeSeed(); + randomizer.Reset(); cout << "# Testing MMFF94 Force Field..." << endl; switch(choice) { diff --git a/tools/obrotamer.cpp b/tools/obrotamer.cpp index 6eb26e35c..133111b82 100644 --- a/tools/obrotamer.cpp +++ b/tools/obrotamer.cpp @@ -75,8 +75,12 @@ int main(int argc,char *argv[]) OBRotorList rl; OBRotamerList rotamers; +#if !OB_USE_OBRANDOMMT OBRandom rand; - rand.TimeSeed(); + rand.Reset(); +#else + OBRandomMT rand{}; +#endif while(ifs.peek() != EOF && ifs.good()) { @@ -95,8 +99,13 @@ int main(int argc,char *argv[]) // (the different angles to sample from the OBRotorRules database) OBRotorIterator ri; OBRotor *rotor = rl.BeginRotor(ri); - for (unsigned int i = 1; i < rl.Size() + 1; ++i, rotor = rl.NextRotor(ri)) - rotorKey[i] = rand.NextInt() % rotor->GetResolution().size(); + for (unsigned int i = 1; i < rl.Size() + 1; ++i, rotor = rl.NextRotor(ri)) { +#if !OB_USE_OBRANDOMMT + rotorKey[i] = rand.UniformInt(0, rotor->GetResolution().size() - 1u); +#else + rotorKey[i] = rand.UniformInt(0, rotor->GetResolution().size() - 1u); +#endif + } rotamers.SetBaseCoordinateSets(mol); rotamers.Setup(mol, rl);