From f5353344b0dd7b6c083340973119282a8b9e2f44 Mon Sep 17 00:00:00 2001 From: Eisuke Kawashima Date: Tue, 26 May 2020 07:44:34 +0900 Subject: [PATCH] Follow the changes of OpenBabel::OBRandom and refine operations related to random numbers --- include/openbabel/conformersearch.h | 9 ++++- include/openbabel/distgeom.h | 6 +++ include/openbabel/forcefield.h | 7 ++++ include/openbabel/math/vector3.h | 11 ++++-- src/conformersearch.cpp | 59 +++++++++++++---------------- src/distgeom.cpp | 15 ++++---- src/forcefield.cpp | 37 +++++++++--------- src/math/vector3.cpp | 30 +++++++-------- test/math.cpp | 10 ++--- tools/obrotamer.cpp | 5 +-- 10 files changed, 103 insertions(+), 86 deletions(-) diff --git a/include/openbabel/conformersearch.h b/include/openbabel/conformersearch.h index ca0f03b208..4165bbb137 100644 --- a/include/openbabel/conformersearch.h +++ b/include/openbabel/conformersearch.h @@ -24,6 +24,7 @@ GNU General Public License for more details. #include #include #include +#include namespace OpenBabel { @@ -31,6 +32,8 @@ namespace OpenBabel { typedef std::vector RotorKeys; typedef std::map,double> mapRotorEnergy; + class OBRandom; + ///@addtogroup conformer Conformer Searching ///@{ @@ -285,6 +288,10 @@ namespace OpenBabel { */ bool Setup(const OBMol &mol, int numConformers = 30, int numChildren = 5, int mutability = 5, int convergence = 25); + /** + * reset prng by the specified seed + */ + void Seed(uint_fast64_t seed); /** * Set the number of conformers. */ @@ -444,7 +451,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 prng; 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/distgeom.h b/include/openbabel/distgeom.h index 93c09c1c1f..4eebbd9c12 100644 --- a/include/openbabel/distgeom.h +++ b/include/openbabel/distgeom.h @@ -24,6 +24,7 @@ GNU General Public License for more details. #include #include +#include #ifndef OBAPI #define OBAPI @@ -38,6 +39,7 @@ namespace OpenBabel { class DistanceGeometryPrivate; class OBCisTransStereo; + class OBRandom; class TetrahedralInfo { int c; @@ -80,6 +82,9 @@ namespace OpenBabel { */ bool Setup(const OBMol &mol, bool useCurrentGeom = false); + //! reset prng by the specified seed + void Seed(uint_fast64_t seed); + void Generate(); void AddConformer(); void GetConformers(OBMol &mol); @@ -115,6 +120,7 @@ namespace OpenBabel { DistanceGeometryPrivate *_d; //!< Internal private data, including bounds matrix Eigen::VectorXd _coord; // one-dimensional vector containing coordinates of atoms std::string input_smiles; + std::unique_ptr prng; unsigned int dim; diff --git a/include/openbabel/forcefield.h b/include/openbabel/forcefield.h index a9fdaa1082..499115074a 100644 --- a/include/openbabel/forcefield.h +++ b/include/openbabel/forcefield.h @@ -19,6 +19,7 @@ GNU General Public License for more details. #ifndef OB_FORCEFIELD_H #define OB_FORCEFIELD_H +#include #include #include @@ -564,6 +565,8 @@ const double GAS_CONSTANT = 8.31446261815324e-3 / KCAL_TO_KJ; //!< kcal mol^-1 std::vector _interGroup; //!< groups for which intra-molecular interactions should be calculated std::vector > _interGroups; //!< groups for which intra-molecular //!< interactions should be calculated + std::shared_ptr prng; + public: /*! Clone the current instance. May be desirable in multithreaded environments, * Should be deleted after use @@ -632,6 +635,10 @@ const double GAS_CONSTANT = 8.31446261815324e-3 / KCAL_TO_KJ; //!< kcal mol^-1 * \return True if successful. */ bool Setup(OBMol &mol, OBFFConstraints &constraints); + /*! + * reset prng by the specified seed + */ + void Seed(uint_fast64_t seed); /*! Load the parameters (this function is overloaded by the individual forcefields, * and is called autoamically from OBForceField::Setup()). */ diff --git a/include/openbabel/math/vector3.h b/include/openbabel/math/vector3.h index 7f8bc5959c..5e1fe8cd58 100644 --- a/include/openbabel/math/vector3.h +++ b/include/openbabel/math/vector3.h @@ -21,9 +21,9 @@ GNU General Public License for more details. #ifndef OB_VECTOR_H #define OB_VECTOR_H -#include -#include +#include #include +#include #include @@ -39,12 +39,14 @@ namespace OpenBabel { class matrix3x3; // declared in math/matrix3x3.h + class OBRandom; // class introduction in vector3.cpp class OBAPI vector3 { private : double _vx, _vy, _vz ; + std::shared_ptr prng; public : //! Constructor @@ -211,6 +213,9 @@ namespace OpenBabel //! \todo Currently unimplemented vector3& operator*= ( const matrix3x3 &); + //! Reset prng by the specified seed + void seed(uint_fast64_t seed); + //! Create a random unit vector void randomUnitVector(); @@ -231,7 +236,7 @@ namespace OpenBabel //! \return The vector length double length () const { - return sqrt( length_2() ); + return std::sqrt( length_2() ); }; //! Access function to get the x-coordinate of the vector const double & x () const diff --git a/src/conformersearch.cpp b/src/conformersearch.cpp index da354f3d3f..8281087de0 100644 --- a/src/conformersearch.cpp +++ b/src/conformersearch.cpp @@ -310,11 +310,7 @@ 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(); + prng.reset(new OBRandom{}); 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 @@ -324,7 +320,6 @@ namespace OpenBabel { OBConformerSearch::~OBConformerSearch() { delete m_filter; - delete (OBRandom*)d; } @@ -390,9 +385,6 @@ namespace OpenBabel { } // create initial population - OBRandom generator; - generator.TimeSeed(); - RotorKey rotorKey(m_rotorList.Size() + 1, 0); // indexed from 1 if (IsGood(rotorKey)) m_rotorKeys.push_back(rotorKey); @@ -408,8 +400,8 @@ 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 (prng->Bernoulli(1.0 / m_mutability)) + rotorKey[i] = prng->UniformInt(0, rotor->GetResolution().size() - 1u); } // duplicates are always rejected if (!IsUniqueKey(m_rotorKeys, rotorKey)) @@ -453,12 +445,14 @@ namespace OpenBabel { return true; } + void OBConformerSearch::Seed(uint_fast64_t seed) + { + prng->Seed(seed); + } + void OBConformerSearch::NextGeneration() { // create next generation population - OBRandom generator; - generator.TimeSeed(); - // generate the children int numConformers = m_rotorKeys.size(); for (int c = 0; c < numConformers; ++c) { @@ -474,8 +468,8 @@ 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 (prng->Bernoulli(1.0 / m_mutability)) + rotorKey[i] = prng->UniformInt(0, rotor->GetResolution().size() - 1u); // permutate gene } // duplicates are always rejected if (!IsUniqueKey(m_rotorKeys, rotorKey)) @@ -789,9 +783,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(); - while (new_val == best[i]) - new_val = ((OBRandom*)d)->NextInt() % rotor->GetResolution().size(); + do { + new_val = prng->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); @@ -835,7 +829,6 @@ namespace OpenBabel { unsigned int i = 0, iniche = 0, j = 0, pop_size = vscores.size (); unsigned int rnd1 = 0, rnd2 = 0, parent1 = 0, parent2 = 0, nsize = 0; int ret_code = 0; - bool flag_crossover = false; OBRotorIterator ri; OBRotor *rotor = nullptr; @@ -843,30 +836,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 = prng->UniformInt(0, pop_size - 1u); + j = prng->UniformInt(0, 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) + const bool flag_crossover = prng->Bernoulli(p_crossover); + if (flag_crossover && prng->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 = prng->UniformInt(0, nsize - 1u); i = dynamic_niches[iniche][rnd1]; - rnd2 = ((OBRandom*)d)->NextInt() % nsize; + rnd2 = prng->UniformInt(0, 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 = prng->UniformInt(0, pop_size - 1u); + j = prng->UniformInt(0, pop_size - 1u); parent2 = vshared_fitnes[i] > vshared_fitnes[j] ? i : j; } @@ -875,7 +868,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 (prng->Bernoulli()) { // Copy parent1 to offspring 1 key1[i] = m_rotorKeys[parent1][i]; key2[i] = m_rotorKeys[parent2][i]; @@ -897,10 +890,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 (prng->Bernoulli(1.0 / m_mutability)) + key1[i] = prng->UniformInt(0, rotor->GetResolution().size() - 1u); + if (prng->Bernoulli(1.0 / m_mutability)) + key2[i] = prng->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 8758ae9e6d..b815d39901 100644 --- a/src/distgeom.cpp +++ b/src/distgeom.cpp @@ -110,7 +110,9 @@ namespace OpenBabel { }; - OBDistanceGeometry::OBDistanceGeometry(): _d(nullptr) {} + OBDistanceGeometry::OBDistanceGeometry(): _d(nullptr) { + prng.reset(new OBRandom{}); + } OBDistanceGeometry::OBDistanceGeometry(const OBMol &mol, bool useCurrentGeometry): _d(nullptr) { @@ -233,6 +235,10 @@ namespace OpenBabel { return true; } + void OBDistanceGeometry::Seed(uint_fast64_t seed) { + prng->Seed(seed); + } + // Set the default bounds to a maximum distance void OBDistanceGeometry::SetUpperBounds() { @@ -1017,13 +1023,11 @@ namespace OpenBabel { unsigned int N = _mol.NumAtoms(); // random distance matrix Eigen::MatrixXd distMat = Eigen::MatrixXd::Zero(N, N); - OBRandom generator; - generator.TimeSeed(); for (size_t i=0; iGetLowerBounds(i, j); double ub = _d->GetUpperBounds(i, j); - double v = generator.NextFloat() * (ub - lb) + lb; + double v = prng->UniformReal(lb, ub); distMat(i, j) = v; distMat(j, i) = v; } @@ -1176,9 +1180,6 @@ namespace OpenBabel { _mol.AddConformer(confCoord); _mol.SetConformer(_mol.NumConformers()); - OBRandom generator(true); // Use system rand() functions - generator.TimeSeed(); - if (_d->debug) { cerr << " max box size: " << _d->maxBoxSize << endl; } diff --git a/src/forcefield.cpp b/src/forcefield.cpp index c8961351af..3f6b576752 100644 --- a/src/forcefield.cpp +++ b/src/forcefield.cpp @@ -931,6 +931,11 @@ namespace OpenBabel return true; } + void OBForceField::Seed(uint_fast64_t seed) + { + prng.reset(new OBRandom{seed}); + } + bool OBForceField::SetLogLevel(int level) { _loglvl = level; @@ -1437,8 +1442,8 @@ namespace OpenBabel OBRotorIterator ri; OBRotor *rotor; - OBRandom generator; - generator.TimeSeed(); + if (!prng) + prng.reset(new OBRandom{}); _origLogLevel = _loglvl; if (_mol.GetCoordinates() == nullptr) @@ -1483,7 +1488,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] = prng->UniformInt(0, rotor->GetResolution().size() - 1u); } rotamers.AddRotamer(rotorKey); } @@ -1607,8 +1612,8 @@ namespace OpenBabel OBRotorIterator ri; OBRotor *rotor; - OBRandom generator; - generator.TimeSeed(); + if (!prng) + prng.reset(new OBRandom{}); int origLogLevel = _loglvl; if (_mol.GetCoordinates() == nullptr) @@ -1732,7 +1737,6 @@ namespace OpenBabel int best_conformer=-1; // double penalty; // for poor performance - double randFloat; // generated random number -- used to pick a rotor double total; // used to calculate the total probability // Start with the current coordinates @@ -1756,19 +1760,18 @@ 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 (prng->Bernoulli(defaultRotor)) // should we just leave this rotor with default setting? continue; - randFloat = generator.NextFloat(); + double p = prng->UniformReal(); total = 0.0; for (unsigned int j = 0; j < rotor->GetResolution().size(); j++) { - if (randFloat > total && randFloat < (total+ rotorWeights[i][j])) { + if (p > total && p < total + rotorWeights[i][j]) { rotorKey[i] = j; break; - } - else + } else { total += rotorWeights[i][j]; + } } } @@ -3421,8 +3424,8 @@ namespace OpenBabel void OBForceField::GenerateVelocities() { cout << "OBForceField::GenerateVelocities()" << endl; - OBRandom generator; - generator.TimeSeed(); + if (!prng) + prng.reset(new OBRandom{}); _ncoords = _mol.NumAtoms() * 3; int velocityIdx; double velocity; @@ -3439,7 +3442,7 @@ namespace OpenBabel if (!_constraints.IsXFixed(a->GetIdx())) { velocity = 0.0; for (int i=0; i < 12; ++i) - velocity += generator.NextFloat(); + velocity += prng->UniformReal(); velocity -= 6.0; velocity *= sqrt((GAS_CONSTANT * _temp)/ (1000 * a->GetAtomicMass())); _velocityPtr[velocityIdx] = velocity; // x10: gromacs uses nm instead of A @@ -3448,7 +3451,7 @@ namespace OpenBabel if (!_constraints.IsYFixed(a->GetIdx())) { velocity = 0.0; for (int i=0; i < 12; ++i) - velocity += generator.NextFloat(); + velocity += prng->UniformReal(); velocity -= 6.0; velocity *= sqrt((GAS_CONSTANT * _temp)/ (1000 * a->GetAtomicMass())); _velocityPtr[velocityIdx+1] = velocity; // idem @@ -3457,7 +3460,7 @@ namespace OpenBabel if (!_constraints.IsZFixed(a->GetIdx())) { velocity = 0.0; for (int i=0; i < 12; ++i) - velocity += generator.NextFloat(); + velocity += prng->UniformReal(); velocity -= 6.0; velocity *= sqrt((GAS_CONSTANT * _temp)/ (1000 * a->GetAtomicMass())); _velocityPtr[velocityIdx+2] = velocity; // idem diff --git a/src/math/vector3.cpp b/src/math/vector3.cpp index e9ca70d1f7..b044af9dd9 100644 --- a/src/math/vector3.cpp +++ b/src/math/vector3.cpp @@ -83,27 +83,23 @@ namespace OpenBabel else return _vz; } - /*! 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. + void vector3::seed(uint_fast64_t seed) + { + prng.reset(new OBRandom{seed}); + } + /*! Replaces *this with a random unit vector, which is (supposed + 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); - l = length_2(); - } - while ( (l > 1.0) || (l < 1e-4) ); - this->normalize(); + if (!prng) + prng.reset(new OBRandom{}); + + double z = prng->UniformReal(-1.0, 1.0); + double phi = prng->UniformReal(0.0, 2 * M_PI); + double rho = sqrt(1.0 - z * z); + this->Set(rho * cos(phi), rho * sin(phi), z); } OBAPI ostream& operator<< ( ostream& co, const vector3& v ) diff --git a/test/math.cpp b/test/math.cpp index a4e18500ab..afd0b8e345 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()); + Diagonal.Set(1, 1, Diagonal.Get(0,0) + randomizer.UniformReal()); + Diagonal.Set(2, 2, Diagonal.Get(1,1) + randomizer.UniformReal()); // 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 b769e0d5fc..cf92d28b21 100644 --- a/tools/obrotamer.cpp +++ b/tools/obrotamer.cpp @@ -75,8 +75,7 @@ int main(int argc,char *argv[]) OBRotorList rl; OBRotamerList rotamers; - OBRandom rand; - rand.TimeSeed(); + OBRandom rand{}; while(ifs.peek() != EOF && ifs.good()) { @@ -96,7 +95,7 @@ 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)) - rotorKey[i] = rand.NextInt() % rotor->GetResolution().size(); + rotorKey[i] = rand.UniformInt(0, rotor->GetResolution().size() - 1u); rotamers.SetBaseCoordinateSets(mol); rotamers.Setup(mol, rl);