From c55eed4b82a8ed4cd414f51df66a61ce27fa10bc Mon Sep 17 00:00:00 2001 From: Eisuke Kawashima Date: Tue, 16 Aug 2022 05:40:49 +0900 Subject: [PATCH 1/9] chore: remove unnecessary forward declaration --- include/openbabel/math/matrix3x3.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/include/openbabel/math/matrix3x3.h b/include/openbabel/math/matrix3x3.h index f50bbc3b90..59e50aa082 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 { From a04c6282304d8ac45fb341210e8c5099ec5c6ab0 Mon Sep 17 00:00:00 2001 From: Eisuke Kawashima Date: Wed, 17 Aug 2022 02:34:18 +0900 Subject: [PATCH 2/9] chore: do not use pointer unnecessarily --- src/math/vector3.cpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/math/vector3.cpp b/src/math/vector3.cpp index a0f79c7622..34e4a3a6bc 100644 --- a/src/math/vector3.cpp +++ b/src/math/vector3.cpp @@ -84,22 +84,18 @@ 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; static OBRandom singleRand(true); - ptr = &singleRand; // 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.NextFloat()-0.5, singleRand.NextFloat()-0.5, singleRand.NextFloat()-0.5); l = length_2(); } while ( (l > 1.0) || (l < 1e-4) ); From f01280e4d20d955c5c7086ac3feab1314da6b9aa Mon Sep 17 00:00:00 2001 From: Eisuke Kawashima Date: Tue, 26 May 2020 07:44:34 +0900 Subject: [PATCH 3/9] refactor: use std::unique_ptr for OpenBabel::OBConformerSearch::d --- include/openbabel/conformersearch.h | 4 ++- src/conformersearch.cpp | 38 +++++++++++++---------------- 2 files changed, 20 insertions(+), 22 deletions(-) diff --git a/include/openbabel/conformersearch.h b/include/openbabel/conformersearch.h index 65bb7ca096..a55ad96705 100644 --- a/include/openbabel/conformersearch.h +++ b/include/openbabel/conformersearch.h @@ -24,8 +24,10 @@ GNU General Public License for more details. #include #include #include +#include namespace OpenBabel { + class OBRandom; typedef std::vector RotorKey; typedef std::vector RotorKeys; @@ -443,7 +445,7 @@ 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 + std::unique_ptr d; 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/src/conformersearch.cpp b/src/conformersearch.cpp index bd26ac0115..53e627993c 100644 --- a/src/conformersearch.cpp +++ b/src/conformersearch.cpp @@ -310,11 +310,8 @@ 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(); + d.reset(new OBRandom()); + d->TimeSeed(); 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 +320,6 @@ namespace OpenBabel { OBConformerSearch::~OBConformerSearch() { - delete (OBRandom*)d; } @@ -792,9 +788,9 @@ 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(); + new_val = d->NextInt() % rotor->GetResolution().size(); while (new_val == best[i]) - new_val = ((OBRandom*)d)->NextInt() % rotor->GetResolution().size(); + new_val = d->NextInt() % rotor->GetResolution().size(); neighbor[i] = new_val; if (IsUniqueKey(backup_population, neighbor) && IsGood(neighbor)) m_rotorKeys.push_back (neighbor); @@ -846,30 +842,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->NextInt() % pop_size; + j = d->NextInt() % pop_size; 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->NextFloat () <= p_crossover); + if (flag_crossover && (d->NextFloat () <= 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->NextInt() % nsize; i = dynamic_niches[iniche][rnd1]; - rnd2 = ((OBRandom*)d)->NextInt() % nsize; + rnd2 = d->NextInt() % nsize; 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->NextInt() % pop_size; + j = d->NextInt() % pop_size; parent2 = vshared_fitnes[i] > vshared_fitnes[j] ? i : j; } @@ -878,7 +874,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->NextInt() % 2) { // Copy parent1 to offspring 1 key1[i] = m_rotorKeys[parent1][i]; key2[i] = m_rotorKeys[parent2][i]; @@ -900,10 +896,10 @@ 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->NextInt() % m_mutability == 0) + key1[i] = d->NextInt() % rotor->GetResolution().size(); + if (d->NextInt() % m_mutability == 0) + key2[i] = d->NextInt() % rotor->GetResolution().size(); } if (IsUniqueKey(m_rotorKeys, key1) && IsGood(key1)) ret_code += 1; From 82c632bb4ff61710fff93fa8c6b8c92d91c44ffd Mon Sep 17 00:00:00 2001 From: Eisuke Kawashima Date: Tue, 16 Aug 2022 09:05:36 +0900 Subject: [PATCH 4/9] feat: add PRNG-related methods to OpenBabel::OBRandom - UniformInt - UniformReal - Normal - Bernoulli --- src/rand.cpp | 13 +++++++++++++ src/rand.h | 14 ++++++++++++++ 2 files changed, 27 insertions(+) diff --git a/src/rand.cpp b/src/rand.cpp index 0975cbf0fb..4630cfb561 100644 --- a/src/rand.cpp +++ b/src/rand.cpp @@ -520,6 +520,19 @@ 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; + } + } //end namespace OpenBabel //! \file rand.cpp diff --git a/src/rand.h b/src/rand.h index 2cb508ed6d..28b48b3c6b 100644 --- a/src/rand.h +++ b/src/rand.h @@ -83,10 +83,24 @@ 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); }; } // end namespace OpenBabel From 29ca772bd456e76c8ae6074ead423f8dbd25a7d5 Mon Sep 17 00:00:00 2001 From: Eisuke Kawashima Date: Tue, 16 Aug 2022 11:01:03 +0900 Subject: [PATCH 5/9] refactor: use newly-added OpenBabel::OBRandom methods --- src/conformersearch.cpp | 50 ++++++++++++++++++++++------------------- src/distgeom.cpp | 6 ++--- src/forcefield.cpp | 36 +++++++++-------------------- src/math/vector3.cpp | 2 +- test/math.cpp | 10 ++++----- tools/obrotamer.cpp | 7 +++--- 6 files changed, 51 insertions(+), 60 deletions(-) diff --git a/src/conformersearch.cpp b/src/conformersearch.cpp index 53e627993c..fdf1556693 100644 --- a/src/conformersearch.cpp +++ b/src/conformersearch.cpp @@ -311,7 +311,7 @@ namespace OpenBabel { niche_mating = 0.7; local_opt_rate = 3; d.reset(new OBRandom()); - d->TimeSeed(); + d->Reset(); 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 @@ -386,7 +386,7 @@ namespace OpenBabel { // create initial population OBRandom generator; - generator.TimeSeed(); + generator.Reset(); RotorKey rotorKey(m_rotorList.Size() + 1, 0); // indexed from 1 if (IsGood(rotorKey)) @@ -403,8 +403,9 @@ 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) { + rotorKey[i] = generator.UniformInt(0, rotor->GetResolution().size() - 1u); + } } // duplicates are always rejected if (!IsUniqueKey(m_rotorKeys, rotorKey)) @@ -452,7 +453,7 @@ namespace OpenBabel { { // create next generation population OBRandom generator; - generator.TimeSeed(); + generator.Reset(); // generate the children int numConformers = m_rotorKeys.size(); @@ -469,8 +470,9 @@ 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) { + rotorKey[i] = generator.UniformInt(0, rotor->GetResolution().size() - 1u); // permutate gene + } } // duplicates are always rejected if (!IsUniqueKey(m_rotorKeys, rotorKey)) @@ -788,9 +790,9 @@ namespace OpenBabel { for (i = 1; i <= m_rotorList.Size(); ++i, rotor = m_rotorList.NextRotor(ri)) { neighbor = best; - new_val = d->NextInt() % rotor->GetResolution().size(); - while (new_val == best[i]) - new_val = d->NextInt() % rotor->GetResolution().size(); + do { + new_val = d->UniformInt(0, rotor->GetResolution().size() - 1u); + } while (new_val == best[i]); neighbor[i] = new_val; if (IsUniqueKey(backup_population, neighbor) && IsGood(neighbor)) m_rotorKeys.push_back (neighbor); @@ -842,30 +844,30 @@ namespace OpenBabel { return 0; // Make a 2-tournament selection to choose first parent - i = d->NextInt() % pop_size; - j = 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 = (d->NextFloat () <= p_crossover); - if (flag_crossover && (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 = d->NextInt() % nsize; + rnd1 = d->UniformInt(0u, nsize - 1u); i = dynamic_niches[iniche][rnd1]; - rnd2 = 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 = d->NextInt() % pop_size; - j = 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; } @@ -874,7 +876,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 (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]; @@ -896,10 +898,12 @@ namespace OpenBabel { rotor = m_rotorList.BeginRotor(ri); for (i = 1; i <= m_rotorList.Size(); ++i, rotor = m_rotorList.NextRotor(ri)) { - if (d->NextInt() % m_mutability == 0) - key1[i] = d->NextInt() % rotor->GetResolution().size(); - if (d->NextInt() % m_mutability == 0) - key2[i] = d->NextInt() % rotor->GetResolution().size(); + if (d->UniformInt(0, m_mutability - 1) == 0) { + key1[i] = d->UniformInt(0, rotor->GetResolution().size() - 1u); + } + if (d->UniformInt(0, m_mutability - 1) == 0) { + key2[i] = d->UniformInt(0, rotor->GetResolution().size() - 1u); + } } if (IsUniqueKey(m_rotorKeys, key1) && IsGood(key1)) ret_code += 1; diff --git a/src/distgeom.cpp b/src/distgeom.cpp index 81fce1c192..fa0fa34f64 100644 --- a/src/distgeom.cpp +++ b/src/distgeom.cpp @@ -1016,12 +1016,12 @@ namespace OpenBabel { // random distance matrix Eigen::MatrixXd distMat = Eigen::MatrixXd::Zero(N, N); OBRandom generator; - generator.TimeSeed(); + generator.Reset(); 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; } @@ -1175,7 +1175,7 @@ namespace OpenBabel { _mol.SetConformer(_mol.NumConformers()); OBRandom generator(true); // Use system rand() functions - generator.TimeSeed(); + generator.Reset(); if (_d->debug) { cerr << " max box size: " << _d->maxBoxSize << endl; diff --git a/src/forcefield.cpp b/src/forcefield.cpp index e314f002a4..0955bd93e6 100644 --- a/src/forcefield.cpp +++ b/src/forcefield.cpp @@ -1434,7 +1434,7 @@ namespace OpenBabel OBRotor *rotor; OBRandom generator; - generator.TimeSeed(); + generator.Reset(); _origLogLevel = _loglvl; if (_mol.GetCoordinates() == nullptr) @@ -1479,7 +1479,7 @@ 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(); + rotorKey[i] = generator.UniformInt(0, rotor->GetResolution().size() - 1u); } rotamers.AddRotamer(rotorKey); } @@ -1604,7 +1604,7 @@ namespace OpenBabel OBRotor *rotor; OBRandom generator; - generator.TimeSeed(); + generator.Reset(); int origLogLevel = _loglvl; if (_mol.GetCoordinates() == nullptr) @@ -1752,11 +1752,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])) { @@ -3417,10 +3416,9 @@ namespace OpenBabel { cout << "OBForceField::GenerateVelocities()" << endl; OBRandom generator; - generator.TimeSeed(); + generator.Reset(); _ncoords = _mol.NumAtoms() * 3; int velocityIdx; - double velocity; _velocityPtr = new double[_ncoords]; memset(_velocityPtr, '\0', sizeof(double)*_ncoords); @@ -3429,32 +3427,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 34e4a3a6bc..afe5b3c4d7 100644 --- a/src/math/vector3.cpp +++ b/src/math/vector3.cpp @@ -95,7 +95,7 @@ namespace OpenBabel double l; do { - this->Set(singleRand.NextFloat()-0.5, singleRand.NextFloat()-0.5, singleRand.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) ); diff --git a/test/math.cpp b/test/math.cpp index 4d776431f1..f0bc26439a 100644 --- a/test/math.cpp +++ b/test/math.cpp @@ -72,7 +72,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 +323,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 +397,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 6eb26e35c5..3a37eb401a 100644 --- a/tools/obrotamer.cpp +++ b/tools/obrotamer.cpp @@ -76,7 +76,7 @@ int main(int argc,char *argv[]) OBRotorList rl; OBRotamerList rotamers; OBRandom rand; - rand.TimeSeed(); + rand.Reset(); while(ifs.peek() != EOF && ifs.good()) { @@ -95,8 +95,9 @@ 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)) { + rotorKey[i] = rand.UniformInt(0, rotor->GetResolution().size() - 1u); + } rotamers.SetBaseCoordinateSets(mol); rotamers.Setup(mol, rl); From d3991533f41d53c097a91cf92cc4fc4e390303ee Mon Sep 17 00:00:00 2001 From: Eisuke Kawashima Date: Mon, 25 May 2020 15:28:38 +0900 Subject: [PATCH 6/9] feat: implement OpenBabel::OBRandomMT as a wrapper of std::mt19937_64 deprecates OpenBabel::OBRandom --- src/rand.cpp | 34 ++++++++++++++++++++++++++++++ src/rand.h | 58 +++++++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 91 insertions(+), 1 deletion(-) diff --git a/src/rand.cpp b/src/rand.cpp index 4630cfb561..5ce264c9ae 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" @@ -533,6 +534,39 @@ namespace OpenBabel 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 28b48b3c6b..9c0f83c3bd 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; @@ -103,6 +105,60 @@ namespace OpenBabel 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 #endif // RAND_H From 4c249105f6bdb89265daf4b7b7fc2c8d8bfe0907 Mon Sep 17 00:00:00 2001 From: Eisuke Kawashima Date: Tue, 16 Aug 2022 12:09:15 +0900 Subject: [PATCH 7/9] refactor: improve method to generate 3-D random unit vector this is enabled if OB_USE_IMPROVED_RANDOM_UNIT_VECTOR, which defaults to whether OpenBabel 4 is released or not, is defined as truthy --- src/math/vector3.cpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/math/vector3.cpp b/src/math/vector3.cpp index afe5b3c4d7..2d5e5d28bf 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 @@ -90,6 +94,7 @@ namespace OpenBabel { static OBRandom singleRand(true); +#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; @@ -100,6 +105,12 @@ namespace OpenBabel } 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 ) From 11629816a51652d9640f0c89e902288a014fc4a8 Mon Sep 17 00:00:00 2001 From: Eisuke Kawashima Date: Tue, 26 May 2020 07:44:34 +0900 Subject: [PATCH 8/9] refactor: switch to OBRandomMT if OB_USE_OBRANDOMMT is defined as truthy the macro defaults to whether OpenBabel 4 is released or not --- include/openbabel/conformersearch.h | 8 ++++++++ src/config.h.cmake | 5 +++++ src/conformersearch.cpp | 32 +++++++++++++++++++++++++++++ src/distgeom.cpp | 8 ++++++++ src/forcefield.cpp | 16 +++++++++++++++ src/math/vector3.cpp | 4 ++++ test/math.cpp | 4 ++++ tools/obrotamer.cpp | 8 ++++++++ 8 files changed, 85 insertions(+) diff --git a/include/openbabel/conformersearch.h b/include/openbabel/conformersearch.h index a55ad96705..733a8e3825 100644 --- a/include/openbabel/conformersearch.h +++ b/include/openbabel/conformersearch.h @@ -27,7 +27,11 @@ GNU General Public License for more details. #include namespace OpenBabel { +#if !OB_USE_OBRANDOMMT class OBRandom; +#else + class OBRandomMT; +#endif typedef std::vector RotorKey; typedef std::vector RotorKeys; @@ -445,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 +#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/src/config.h.cmake b/src/config.h.cmake index b2977e4ccc..b438dab1d2 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 fdf1556693..c917b73f38 100644 --- a/src/conformersearch.cpp +++ b/src/conformersearch.cpp @@ -310,8 +310,12 @@ namespace OpenBabel { p_crossover = 0.7; niche_mating = 0.7; local_opt_rate = 3; +#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 @@ -385,8 +389,12 @@ namespace OpenBabel { } // create initial population +#if !OB_USE_OBRANDOMMT OBRandom generator; generator.Reset(); +#else + OBRandomMT generator{}; +#endif RotorKey rotorKey(m_rotorList.Size() + 1, 0); // indexed from 1 if (IsGood(rotorKey)) @@ -404,7 +412,11 @@ namespace OpenBabel { OBRotor *rotor = m_rotorList.BeginRotor(ri); for (unsigned int i = 1; i < m_rotorList.Size() + 1; ++i, rotor = m_rotorList.NextRotor(ri)) { 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 @@ -452,8 +464,12 @@ namespace OpenBabel { void OBConformerSearch::NextGeneration() { // create next generation population +#if !OB_USE_OBRANDOMMT OBRandom generator; generator.Reset(); +#else + OBRandomMT generator{}; +#endif // generate the children int numConformers = m_rotorKeys.size(); @@ -471,7 +487,11 @@ namespace OpenBabel { OBRotor *rotor = m_rotorList.BeginRotor(ri); for (unsigned int i = 1; i < m_rotorList.Size() + 1; ++i, rotor = m_rotorList.NextRotor(ri)) { 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 @@ -791,7 +811,11 @@ namespace OpenBabel { { neighbor = best; 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)) @@ -899,10 +923,18 @@ namespace OpenBabel { for (i = 1; i <= m_rotorList.Size(); ++i, rotor = m_rotorList.NextRotor(ri)) { 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)) diff --git a/src/distgeom.cpp b/src/distgeom.cpp index fa0fa34f64..4c2bd01917 100644 --- a/src/distgeom.cpp +++ b/src/distgeom.cpp @@ -1015,8 +1015,12 @@ 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.Reset(); +#else + OBRandomMT generator{}; +#endif for (size_t i=0; iGetLowerBounds(i, j); @@ -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.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 0955bd93e6..889fd29fb0 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.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 +#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.Reset(); +#else + OBRandomMT generator{}; +#endif int origLogLevel = _loglvl; if (_mol.GetCoordinates() == nullptr) @@ -3415,8 +3427,12 @@ namespace OpenBabel void OBForceField::GenerateVelocities() { cout << "OBForceField::GenerateVelocities()" << endl; +#if !OB_USE_OBRANDOMMT OBRandom generator; generator.Reset(); +#else + OBRandomMT generator{}; +#endif _ncoords = _mol.NumAtoms() * 3; int velocityIdx; diff --git a/src/math/vector3.cpp b/src/math/vector3.cpp index 2d5e5d28bf..3fe1343b2e 100644 --- a/src/math/vector3.cpp +++ b/src/math/vector3.cpp @@ -92,7 +92,11 @@ namespace OpenBabel */ void vector3::randomUnitVector() { +#if !OB_USE_OBRANDOMMT static OBRandom singleRand(true); +#else + static OBRandomMT singleRand{}; +#endif #if !OB_USE_IMPROVED_RANDOM_UNIT_VECTOR // obtain a random vector with 0.001 <= length^2 <= 1.0, normalize diff --git a/test/math.cpp b/test/math.cpp index f0bc26439a..bbe815303a 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]; diff --git a/tools/obrotamer.cpp b/tools/obrotamer.cpp index 3a37eb401a..133111b821 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.Reset(); +#else + OBRandomMT rand{}; +#endif while(ifs.peek() != EOF && ifs.good()) { @@ -96,7 +100,11 @@ int main(int argc,char *argv[]) OBRotorIterator ri; OBRotor *rotor = rl.BeginRotor(ri); 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); From 632eaa7e449c7774cb7ad2b51629ea486c347fcf Mon Sep 17 00:00:00 2001 From: Eisuke Kawashima Date: Wed, 17 Aug 2022 11:10:24 +0900 Subject: [PATCH 9/9] docs: add explanations of OB_RANDOM_SEED environment variable --- doc/obabel.1 | 7 +++++++ doc/obdistgen.1 | 7 +++++++ doc/obgen.1 | 7 +++++++ doc/obrotamer.1 | 7 +++++++ 4 files changed, 28 insertions(+) diff --git a/doc/obabel.1 b/doc/obabel.1 index 1411ccca14..f576bec0e7 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 786acf813a..409b0f5666 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 7901865a51..9ba2d9a4c2 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 adc47e4bad..b8f1182a16 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