ATOOLS/Math/Random.H

0001 #ifndef ATOOLS_Math_Random_H
0002 #define ATOOLS_Math_Random_H
0003
0004 #include <limits>
0005 #include <fstream>
0006 #include <cmath>
0007 #include <algorithm>
0008 #include <stddef.h>
0009 #include "ATOOLS/Math/MathTools.H"
0010 #include "ATOOLS/Org/Terminator_Objects.H"
0011 #include "ATOOLS/Org/Getter_Function.H"
0012
0013 namespace ATOOLS {
0014
0015   struct RNG_Key {
0016   };
0017
0018   class External_RNG {
0019   public:
0020
0021     virtual ~External_RNG();
0022
0023     virtual double Get() = 0;
0024
0025     // Derived classes can opt in to support status saving/restoring.
0026     // In that case, override all three of the below functions and
0027     // in particular return true in CanRestoreStatus().
0028     virtual bool CanRestoreStatus() const { return false; }
0029     virtual void SaveStatus() {}
0030     virtual void RestoreStatus() {}
0031
0032   };
0033
0034   typedef Getter_Function<External_RNG,RNG_Key> RNG_Getter;
0035
0036   class Marsaglia;
0037
0038   class Random: public Terminator_Object {
0039   private:
0040
0041     int      activeGenerator;
0042
0043     long int m_id, m_sid;
0044
0045     std::stringstream m_lastincrementedseed;
0046     size_t m_nsinceinit, m_increment;
0047
0048     External_RNG *p_external;
0049
0050     Marsaglia *p_ran4[2];
0051
0052     double Ran2(long *idum);
0053     double Ran4();
0054
0055     bool ReadInStatus(const std::string &path);
0056     void PrepareTerminate();
0057
0058     // temporary methods for Ran4()
0059     int  WriteOutStatus4(const char *outfile);
0060     int  WriteOutStatus4(std::ostream &os,const size_t &idx);
0061     int  WriteOutSavedStatus4(const char *outfile);
0062     void ReadInStatus4(const char * filename);
0063     size_t ReadInStatus4(std::istream &is,const size_t &idx);
0064     void SaveStatus4();
0065     void RestoreStatus4();
0066
0067   public:
0068
0069     // constructors
0070     Random(long nid);  // initialization for Ran2()
0071     Random(unsigned int i1,unsigned int i2,unsigned int i3,
0072        unsigned int i4,unsigned int i5,unsigned int i6);
0073
0074     // destructor
0075     ~Random();
0076
0077     // member functions
0078     bool InitExternal();
0079     void SetSeed(long nid);  // seed for Rnd2()
0080     void SetSeed(unsigned int i1,unsigned int i2,
0081          unsigned int i3,unsigned int i4);
0082     long int GetSeed() { return m_id; }
0083
0084     int  WriteOutStatus(const char *outfile);
0085     int  WriteOutStatus(std::ostream &os,const size_t &idx);
0086     int  WriteOutSavedStatus(const char *outfile);
0087     void ReadInStatus(const char *infile);
0088     size_t ReadInStatus(std::istream &is,const size_t &idx);
0089     bool CanRestoreStatus() const;
0090     void SaveStatus();
0091     void RestoreStatus();
0092     void FastForward(const size_t &n);
0093     void ResetToLastIncrementedSeed();
0094     void EraseLastIncrementedSeed()
0095     { m_lastincrementedseed.str(std::string()); }
0096     void SetSeedStorageIncrement(size_t inc) { m_increment=inc; }
0097
0098     // return uniformly distributed random number in [0,1] using active Generator
0099     double Get();
0100     // produce Gaussian distributed random number using active Generator
0101     // according to the Marsaglia method
0102     double GetGaussian() {
0103       static auto hasSpare = false;
0104       static double spare;
0105       if (hasSpare) {
0106         hasSpare = false;
0107         return spare;
0108       }
0109       double ran1, ran2, R;
0110       do {
0111         ran1 = 2.*Get()-1.;
0112         ran2 = 2.*Get()-1.;
0113         R    = ran1*ran1+ran2*ran2;
0114       } while (R>1. || R==0.);
0115       R  = std::sqrt(-2.*std::log(R)/R);
0116       hasSpare = true;
0117       spare = ran2 * R;
0118       return ran1 * R;
0119     }
0120     // produce Poissonian distributed random number using active Generator
0121     double Poissonian(const double & lambda) {
0122       if(lambda>500.) {
0123     double u = Get();
0124     double v = Get();
0125     return int(std::sqrt(lambda)*std::sqrt(-2.*std::log(u))*std::cos(2.*M_PI*v)+lambda);
0126       }
0127       double disc(std::exp(-lambda)),p(1.);
0128       int N(0);
0129       while ((p*=Get())>disc) N++;
0130       return N;
0131     }
0132     double Theta() { return std::acos(2.*Get()-1.); }
0133     double GetNZ();
0134
0135     External_RNG* GetExternalRng() { return p_external; }
0136
0137     // adopt requirements for UniformRandomBitGenerator
0138     // (e.g. for use in std::shuffle)
0139     typedef size_t result_type;
0140     result_type operator()();
0141     constexpr static result_type min() noexcept {
0142       return std::numeric_limits<result_type>::min();
0143     }
0144     constexpr static result_type max() noexcept {
0145       return std::numeric_limits<result_type>::max();
0146     }
0147
0148   };// end of class Random
0149
0150   extern Random *ran;
0151
0152   // --------------------------------------------------
0153   //         Doxygen part
0154   // --------------------------------------------------
0155
0156   /*!
0157     \file
0158     \brief contains the class Random
0159   */
0160
0161   /*!
0162     \class Random
0163     \brief supplies uniformly distributed random numbers
0164   */
0165
0166   /*!
0167     \fn double Random::Ran2(long *idum)
0168     \brief is a very stable and powerful random number routine
0169   */
0170
0171   /*!
0172     \fn double Random::Ran4()
0173     \brief a new random generator that still needs to be tested
0174   */
0175
0176   /*!
0177     \fn Random::Random(long nid)
0178     \brief Constructor initialises the random number generator with a given seed
0179   */
0180
0181   /*!
0182     \fn Random::Random(unsigned int i1, unsigned int i2, unsigned int i3, unsigned int i4, unsigned int i5, unsigned int i6)
0183     \brief A constructor that initializes the Rnd4() routine.
0184
0185     Even though there are two different constructors for Rnd2() and Rnd4(),
0186     it is possible to switch between the two routines by calling the
0187     corresponding SetSeed() method.
0188   */
0189
0190   /*!
0191     \fn Random::~Random()
0192     \brief Destructor
0193   */
0194
0195   /*!
0196     \fn double Random::Get()
0197     \brief is the main routine, returns a single random number in [0,1]
0198
0199     The number is determined either by using Ran2() or Ran4(), depending on
0200     which of the two generators is set in the activeGenerator variable.
0201   */
0202
0203   /*!
0204     \fn double Random::GetNZ()
0205     \brief retrun a not zero random number
0206   */
0207
0208   /*!
0209     \fn long Random::GetSeed()
0210     \brief returns a the seed
0211
0212     No corresponding method for Ran4() exists so far.
0213   */
0214
0215   /*!
0216     \fn void Random::SetSeed(long nid)
0217     \brief sets a new seed and (re)initializes the random number generator Rnd2()
0218   */
0219
0220   /*!
0221     \fn void Random::SetSeed(unsigned int i1, unsigned int i2, unsigned int i3, unsigned int i4)
0222     \brief sets a new seed and (re)initializes the random number generator Rnd4()
0223   */
0224
0225   /*!
0226     \fn double Random::Theta()
0227     \brief returns an angle \f$\phi\f$ for a uniform \f$cos(\phi)\f$ distribution
0228   */
0229
0230   /*!
0231     \fn int Random::WriteOutStatus(const char* outfile)
0232     \brief writes the complete status the random generator in a file
0233
0234     This method can be used to save the status of a random generator in a file
0235     the number of its entry in this file is return and can be used to read in
0236     the status via \link ReadInStatus \endlink.
0237   */
0238
0239   /*!
0240     \fn bool Random::ReadInStatus(const char* infile)
0241     \brief reads in a status from a file
0242   */
0243
0244 }// end of namespace ATOOLS
0245
0246 #endif