Skip to content

Commit

Permalink
Follow the changes of OpenBabel::OBRandom and refine operations relat…
Browse files Browse the repository at this point in the history 
…ed to random numbers
  • Loading branch information
@e-kwsm
e-kwsm committed May 26, 2020
1 parent 2965884 commit f535334
Show file tree
Hide file tree
Showing 10 changed files with 103 additions and 86 deletions.
9 changes: 8 additions & 1 deletion include/openbabel/conformersearch.h
View file
Expand Up @@ -24,13 +24,16 @@ GNU General Public License for more details.
#include <openbabel/mol.h>
#include <openbabel/rotor.h>
#include <openbabel/rotamer.h>
#include <memory>

namespace OpenBabel {

typedef std::vector<int> RotorKey;
typedef std::vector<RotorKey> RotorKeys;
typedef std::map<std::vector<int>,double> mapRotorEnergy;

class OBRandom;

///@addtogroup conformer Conformer Searching
///@{

Expand Down Expand Up @@ -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.
*/
Expand Down Expand Up @@ -444,7 +451,7 @@ 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
std::unique_ptr<OBRandom> 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
Expand Down
6 changes: 6 additions & 0 deletions include/openbabel/distgeom.h
View file
Expand Up @@ -24,6 +24,7 @@ GNU General Public License for more details.
#include <openbabel/mol.h>

#include <iostream>
#include <memory>

#ifndef OBAPI
#define OBAPI
Expand All @@ -38,6 +39,7 @@ namespace OpenBabel {

class DistanceGeometryPrivate;
class OBCisTransStereo;
class OBRandom;

class TetrahedralInfo {
int c;
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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<OBRandom> prng;

unsigned int dim;

Expand Down
7 changes: 7 additions & 0 deletions include/openbabel/forcefield.h
View file
Expand Up @@ -19,6 +19,7 @@ GNU General Public License for more details.
#ifndef OB_FORCEFIELD_H
#define OB_FORCEFIELD_H

#include <memory>
#include <vector>
#include <string>

Expand Down Expand Up @@ -564,6 +565,8 @@ const double GAS_CONSTANT = 8.31446261815324e-3 / KCAL_TO_KJ; //!< kcal mol^-1
std::vector<OBBitVec> _interGroup; //!< groups for which intra-molecular interactions should be calculated
std::vector<std::pair<OBBitVec, OBBitVec> > _interGroups; //!< groups for which intra-molecular
//!< interactions should be calculated
std::shared_ptr<OBRandom> prng;

public:
/*! Clone the current instance. May be desirable in multithreaded environments,
* Should be deleted after use
Expand Down Expand Up @@ -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()).
*/
Expand Down
11 changes: 8 additions & 3 deletions include/openbabel/math/vector3.h
View file
Expand Up @@ -21,9 +21,9 @@ GNU General Public License for more details.
#ifndef OB_VECTOR_H
#define OB_VECTOR_H

#include <ostream>
#include <math.h>
#include <cmath>
#include <iostream>
#include <memory>

#include <openbabel/babelconfig.h>

Expand All @@ -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<OBRandom> prng;

public :
//! Constructor
Expand Down Expand Up @@ -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();

Expand All @@ -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
Expand Down
59 changes: 26 additions & 33 deletions src/conformersearch.cpp
View file
Expand Up @@ -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
Expand All @@ -324,7 +320,6 @@ namespace OpenBabel {
OBConformerSearch::~OBConformerSearch()
{
delete m_filter;
delete (OBRandom*)d;
}


Expand Down Expand Up @@ -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);
Expand All @@ -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))
Expand Down Expand Up @@ -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) {
Expand All @@ -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))
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -835,38 +829,37 @@ 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;

if (pop_size < 2)
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;
}

Expand All @@ -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];
Expand All @@ -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;
Expand Down
15 changes: 8 additions & 7 deletions src/distgeom.cpp
View file
Expand Up @@ -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)
{
Expand Down Expand Up @@ -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()
{
Expand Down Expand Up @@ -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; 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 = prng->UniformReal(lb, ub);
distMat(i, j) = v;
distMat(j, i) = v;
}
Expand Down Expand Up @@ -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;
}
Expand Down

0 comments on commit f535334

Please sign in to comment.