Random.h

#ifndef UTIL_RANDOM_H
#define UTIL_RANDOM_H

/*
* Util Package - C++ Utilities for Scientific Computation
*
* Copyright 2010 - 2017, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/

#include <util/param/ParamComposite.h>
#include <util/space/Vector.h>

#include <cmath>

#include <util/random/mersenne/mtrand.h>
#define UTIL_ENGINE MTRand_int32

namespace Util
{

   class Vector;

   class Random : public ParamComposite
   {
 
   public:

      typedef unsigned long SeedType;

      Random();

      virtual ~Random();

      virtual void readParameters(std::istream &in);
   
      virtual void loadParameters(Serializable::IArchive &ar);
   
      virtual void save(Serializable::OArchive &ar);
   
      void setSeed(SeedType seed);
   
      double uniform();
   
      double uniform(double range1, double range2);
   
      long uniformInt(long range1, long range2);
    
      void getPoint(double minR[], double maxR[], double r[]);
   
      double gaussian(void);
   
      void unitVector(Vector& v);
   
      bool metropolis(double ratio);

      long drawFrom(double probability[], long size);
    
      template <class Archive>
      void serialize(Archive& ar, const unsigned int version);

      long seed();
   
   private:

      UTIL_ENGINE engine_;

      SeedType seed_;

      bool isInitialized_;
   
      void setSeed();
   
   };
  
   // Inline methods 

   /*
   * Get a uniform double precision number.
   */
   inline double Random::uniform()
   {
      // Divide 32 bit integer engine output by 2^32
      return static_cast<double>(engine_()) * (1.0/ 4294967296.0); 
   }
 
   /*
   * Get a uniform double precision number.
   */
   inline double Random::uniform(double range1, double range2)
   {
      // Divide 32 bit integer engine output by 2^32
      double frac = static_cast<double>(engine_()) * (1.0/ 4294967296.0); 
      return range1 + (range2 - range1)*frac; 
   }
 
   /* 
   * Return a random long integer x uniformly distributed in range1 <= x < range2.
   *
   * Parameters range1 and range2 must be within the range of long integers.
   */
   inline long Random::uniformInt(long range1, long range2) {
      double x = uniform(range1, range2);
      return long(x);
   }
   
   /* 
   * Choose one of several outcomes with a specified set of probabilities.
   *
   * Precondition: Elements of probability array must add to 1.0
   *
   * \param probability[] array of probabilities, for indices 0,...,size-1
   * \param size          number of options
   */
   inline long Random::drawFrom(double probability[], long size) {
      double roulette = uniform();
      long   n=0;
      double cumProb = probability[0];
      while ( cumProb < roulette && n < size ) {
         ++n;
         cumProb += probability[n];
      }
      return(n);
   }
   
   /*
   * Get random point within a rectangular prism.
   */
   inline void Random::getPoint(double minR[], double maxR[], double r[]) {
      r[0] = uniform(minR[0], maxR[0]);
      r[1] = uniform(minR[1], maxR[1]);
      r[2] = uniform(minR[2], maxR[2]);
   }
   
   /* 
   * Metropolis criteria for whether to accept a MC move.
   */
   inline bool Random::metropolis(double ratio)
   {
      if ( ratio > 1.0 ) {
         return true;
      } else {
         double ran = uniform();
         if ( ran < ratio ) {
            return true;
         } else {
            return false;
         }
      }
   }

   /* 
   * Returns value of random seed (private member variable idum)
   */
   inline long Random::seed() 
   {  return seed_; }

   /*
   * Serialize to/from an archive.
   */
   template <class Archive>
   void Random::serialize(Archive& ar, const unsigned int version)
   {
      ar & engine_;
      ar & seed_;
   }

} 
#undef UTIL_ENGINE 
#endif // ifndef UTIL_RANDOM_H