Reimplement OpenBabel::OBRandom as a wrapper of std::mt19937_64

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"

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

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:

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

include/openbabel/conformersearch.h

@@ -24,8 +24,14 @@ GNU General Public License for more details.
 #include <openbabel/mol.h>
 #include <openbabel/rotor.h>
 #include <openbabel/rotamer.h>
+#include <memory>
 
 namespace OpenBabel {
+#if !OB_USE_OBRANDOMMT
+  class OBRandom;
+#else
+  class OBRandomMT;
+#endif
 
   typedef std::vector<int> RotorKey;
   typedef std::vector<RotorKey> RotorKeys;
@@ -443,7 +449,11 @@ namespace OpenBabel {
       std::vector<std::vector <int> > dynamic_niches; //!< The dynamically found niches
       std::vector<int> 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<OBRandom> d;
+#else
+      std::unique_ptr<OBRandomMT> 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

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
     {

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

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<int>(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<int>(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<int>(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<int>(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<int>(0, rotor->GetResolution().size() - 1u);
+#endif
+        }
       }
     if (IsUniqueKey(m_rotorKeys, key1) && IsGood(key1))
       ret_code += 1;

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; i<N; ++i) {
       for(size_t j=0; j<i; ++j) {
         double lb = _d->GetLowerBounds(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;