COIN-OR::LEMON - Graph Library

Changeset 1404:c8d0179a32a2 in lemon for lemon/random.h


Ignore:
Timestamp:
10/17/18 19:22:52 (6 years ago)
Author:
Alpar Juttner <alpar@…>
Branch:
default
Parents:
1389:57167d92e96c (diff), 1402:3c00344f49c9 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Phase:
public
Message:

Merge bugfixes #610,#611,#612,#614

Files:
2 edited

Legend:

Unmodified
Added
Removed
  • lemon/random.h

    r1380 r1404  
    340340          num = rnd() & mask;
    341341        } while (num > max);
    342         return num;
     342        return static_cast<Result>(num);
    343343      }
    344344    };
  • lemon/random.h

    r1396 r1404  
    111111      static const Word loMask = (1u << 31) - 1;
    112112      static const Word hiMask = ~loMask;
    113 
    114113
    115114      static Word tempering(Word rnd) {
     
    244243    private:
    245244
    246 
    247245      void fillState() {
    248246        static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask };
     
    252250        current = state + length;
    253251
    254         register Word *curr = state + length - 1;
    255         register long num;
     252        Word *curr = state + length - 1;
     253        long num;
    256254
    257255        num = length - shift;
     
    272270      }
    273271
    274 
    275272      Word *current;
    276273      Word state[length];
     
    297294    template <typename Result, typename Word,
    298295              int rest = std::numeric_limits<Result>::digits, int shift = 0,
    299               bool last = rest <= std::numeric_limits<Word>::digits>
     296              bool last = (rest <= std::numeric_limits<Word>::digits)>
    300297    struct IntConversion {
    301298      static const int bits = std::numeric_limits<Word>::digits;
     
    466463    };
    467464
    468   }
     465    /// \ingroup misc
     466    ///
     467    /// \brief Mersenne Twister random number generator
     468    ///
     469    /// The Mersenne Twister is a twisted generalized feedback
     470    /// shift-register generator of Matsumoto and Nishimura. The period
     471    /// of this generator is \f$ 2^{19937} - 1\f$ and it is
     472    /// equi-distributed in 623 dimensions for 32-bit numbers. The time
     473    /// performance of this generator is comparable to the commonly used
     474    /// generators.
     475    ///
     476    /// This is a template implementation of both 32-bit and
     477    /// 64-bit architecture optimized versions. The generators differ
     478    /// sligthly in the initialization and generation phase so they
     479    /// produce two completly different sequences.
     480    ///
     481    /// \alert Do not use this class directly, but instead one of \ref
     482    /// Random, \ref Random32 or \ref Random64.
     483    ///
     484    /// The generator gives back random numbers of serveral types. To
     485    /// get a random number from a range of a floating point type, you
     486    /// can use one form of the \c operator() or the \c real() member
     487    /// function. If you want to get random number from the {0, 1, ...,
     488    /// n-1} integer range, use the \c operator[] or the \c integer()
     489    /// method. And to get random number from the whole range of an
     490    /// integer type, you can use the argumentless \c integer() or
     491    /// \c uinteger() functions. Finally, you can get random bool with
     492    /// equal chance of true and false or with given probability of true
     493    /// result using the \c boolean() member functions.
     494    ///
     495    /// Various non-uniform distributions are also supported: normal (Gauss),
     496    /// exponential, gamma, Poisson, etc.; and a few two-dimensional
     497    /// distributions, too.
     498    ///
     499    ///\code
     500    /// // The commented code is identical to the other
     501    /// double a = rnd();                     // [0.0, 1.0)
     502    /// // double a = rnd.real();             // [0.0, 1.0)
     503    /// double b = rnd(100.0);                // [0.0, 100.0)
     504    /// // double b = rnd.real(100.0);        // [0.0, 100.0)
     505    /// double c = rnd(1.0, 2.0);             // [1.0, 2.0)
     506    /// // double c = rnd.real(1.0, 2.0);     // [1.0, 2.0)
     507    /// int d = rnd[100000];                  // 0..99999
     508    /// // int d = rnd.integer(100000);       // 0..99999
     509    /// int e = rnd[6] + 1;                   // 1..6
     510    /// // int e = rnd.integer(1, 1 + 6);     // 1..6
     511    /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1
     512    /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1
     513    /// bool g = rnd.boolean();               // P(g = true) = 0.5
     514    /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8
     515    ///\endcode
     516    ///
     517    /// LEMON provides a global instance of the random number generator:
     518    /// \ref lemon::rnd "rnd". In most cases, it is a good practice
     519    /// to use this global generator to get random numbers.
     520    ///
     521    /// \sa \ref Random, \ref Random32 or \ref Random64.
     522    template<class Word>
     523    class Random {
     524    private:
     525
     526      _random_bits::RandomCore<Word> core;
     527      _random_bits::BoolProducer<Word> bool_producer;
     528
     529
     530    public:
     531
     532      ///\name Initialization
     533      ///
     534      /// @{
     535
     536      /// \brief Default constructor
     537      ///
     538      /// Constructor with constant seeding.
     539      Random() { core.initState(); }
     540
     541      /// \brief Constructor with seed
     542      ///
     543      /// Constructor with seed. The current number type will be converted
     544      /// to the architecture word type.
     545      template <typename Number>
     546      Random(Number seed) {
     547        _random_bits::Initializer<Number, Word>::init(core, seed);
     548      }
     549
     550      /// \brief Constructor with array seeding
     551      ///
     552      /// Constructor with array seeding. The given range should contain
     553      /// any number type and the numbers will be converted to the
     554      /// architecture word type.
     555      template <typename Iterator>
     556      Random(Iterator begin, Iterator end) {
     557        typedef typename std::iterator_traits<Iterator>::value_type Number;
     558        _random_bits::Initializer<Number, Word>::init(core, begin, end);
     559      }
     560
     561      /// \brief Copy constructor
     562      ///
     563      /// Copy constructor. The generated sequence will be identical to
     564      /// the other sequence. It can be used to save the current state
     565      /// of the generator and later use it to generate the same
     566      /// sequence.
     567      Random(const Random& other) {
     568        core.copyState(other.core);
     569      }
     570
     571      /// \brief Assign operator
     572      ///
     573      /// Assign operator. The generated sequence will be identical to
     574      /// the other sequence. It can be used to save the current state
     575      /// of the generator and later use it to generate the same
     576      /// sequence.
     577      Random& operator=(const Random& other) {
     578        if (&other != this) {
     579          core.copyState(other.core);
     580        }
     581        return *this;
     582      }
     583
     584      /// \brief Seeding random sequence
     585      ///
     586      /// Seeding the random sequence. The current number type will be
     587      /// converted to the architecture word type.
     588      template <typename Number>
     589      void seed(Number seed) {
     590        _random_bits::Initializer<Number, Word>::init(core, seed);
     591      }
     592
     593      /// \brief Seeding random sequence
     594      ///
     595      /// Seeding the random sequence. The given range should contain
     596      /// any number type and the numbers will be converted to the
     597      /// architecture word type.
     598      template <typename Iterator>
     599      void seed(Iterator begin, Iterator end) {
     600        typedef typename std::iterator_traits<Iterator>::value_type Number;
     601        _random_bits::Initializer<Number, Word>::init(core, begin, end);
     602      }
     603
     604      /// \brief Seeding from file or from process id and time
     605      ///
     606      /// By default, this function calls the \c seedFromFile() member
     607      /// function with the <tt>/dev/urandom</tt> file. If it does not success,
     608      /// it uses the \c seedFromTime().
     609      /// \return Currently always \c true.
     610      bool seed() {
     611#ifndef LEMON_WIN32
     612        if (seedFromFile("/dev/urandom", 0)) return true;
     613#endif
     614        if (seedFromTime()) return true;
     615        return false;
     616      }
     617
     618      /// \brief Seeding from file
     619      ///
     620      /// Seeding the random sequence from file. The linux kernel has two
     621      /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which
     622      /// could give good seed values for pseudo random generators (The
     623      /// difference between two devices is that the <tt>random</tt> may
     624      /// block the reading operation while the kernel can give good
     625      /// source of randomness, while the <tt>urandom</tt> does not
     626      /// block the input, but it could give back bytes with worse
     627      /// entropy).
     628      /// \param file The source file
     629      /// \param offset The offset, from the file read.
     630      /// \return \c true when the seeding successes.
     631#ifndef LEMON_WIN32
     632      bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)
     633#else
     634        bool seedFromFile(const std::string& file = "", int offset = 0)
     635#endif
     636      {
     637        std::ifstream rs(file.c_str());
     638        const int size = 4;
     639        Word buf[size];
     640        if (offset != 0 && !rs.seekg(offset)) return false;
     641        if (!rs.read(reinterpret_cast<char*>(buf), sizeof(buf))) return false;
     642        seed(buf, buf + size);
     643        return true;
     644      }
     645
     646      /// \brief Seeding from process id and time
     647      ///
     648      /// Seeding from process id and time. This function uses the
     649      /// current process id and the current time for initialize the
     650      /// random sequence.
     651      /// \return Currently always \c true.
     652      bool seedFromTime() {
     653#ifndef LEMON_WIN32
     654        timeval tv;
     655        gettimeofday(&tv, 0);
     656        seed(getpid() + tv.tv_sec + tv.tv_usec);
     657#else
     658        seed(bits::getWinRndSeed());
     659#endif
     660        return true;
     661      }
     662
     663      /// @}
     664
     665      ///\name Uniform Distributions
     666      ///
     667      /// @{
     668
     669      /// \brief Returns a random real number from the range [0, 1)
     670      ///
     671      /// It returns a random real number from the range [0, 1). The
     672      /// default Number type is \c double.
     673      template <typename Number>
     674      Number real() {
     675        return _random_bits::RealConversion<Number, Word>::convert(core);
     676      }
     677
     678      double real() {
     679        return real<double>();
     680      }
     681
     682      /// \brief Returns a random real number from the range [0, 1)
     683      ///
     684      /// It returns a random double from the range [0, 1).
     685      double operator()() {
     686        return real<double>();
     687      }
     688
     689      /// \brief Returns a random real number from the range [0, b)
     690      ///
     691      /// It returns a random real number from the range [0, b).
     692      double operator()(double b) {
     693        return real<double>() * b;
     694      }
     695
     696      /// \brief Returns a random real number from the range [a, b)
     697      ///
     698      /// It returns a random real number from the range [a, b).
     699      double operator()(double a, double b) {
     700        return real<double>() * (b - a) + a;
     701      }
     702
     703      /// \brief Returns a random integer from a range
     704      ///
     705      /// It returns a random integer from the range {0, 1, ..., b - 1}.
     706      template <typename Number>
     707      Number integer(Number b) {
     708        return _random_bits::Mapping<Number, Word>::map(core, b);
     709      }
     710
     711      /// \brief Returns a random integer from a range
     712      ///
     713      /// It returns a random integer from the range {a, a + 1, ..., b - 1}.
     714      template <typename Number>
     715      Number integer(Number a, Number b) {
     716        return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;
     717      }
     718
     719      /// \brief Returns a random integer from a range
     720      ///
     721      /// It returns a random integer from the range {0, 1, ..., b - 1}.
     722      template <typename Number>
     723      Number operator[](Number b) {
     724        return _random_bits::Mapping<Number, Word>::map(core, b);
     725      }
     726
     727      /// \brief Returns a random non-negative integer
     728      ///
     729      /// It returns a random non-negative integer uniformly from the
     730      /// whole range of the current \c Number type. The default result
     731      /// type of this function is <tt>unsigned int</tt>.
     732      template <typename Number>
     733      Number uinteger() {
     734        return _random_bits::IntConversion<Number, Word>::convert(core);
     735      }
     736
     737      unsigned int uinteger() {
     738        return uinteger<unsigned int>();
     739      }
     740
     741      /// \brief Returns a random integer
     742      ///
     743      /// It returns a random integer uniformly from the whole range of
     744      /// the current \c Number type. The default result type of this
     745      /// function is \c int.
     746      template <typename Number>
     747      Number integer() {
     748        static const int nb = std::numeric_limits<Number>::digits +
     749          (std::numeric_limits<Number>::is_signed ? 1 : 0);
     750        return _random_bits::IntConversion<Number, Word, nb>::convert(core);
     751      }
     752
     753      int integer() {
     754        return integer<int>();
     755      }
     756
     757      /// \brief Returns a random bool
     758      ///
     759      /// It returns a random bool. The generator holds a buffer for
     760      /// random bits. Every time when it become empty the generator makes
     761      /// a new random word and fill the buffer up.
     762      bool boolean() {
     763        return bool_producer.convert(core);
     764      }
     765
     766      /// @}
     767
     768      ///\name Non-uniform Distributions
     769      ///
     770      ///@{
     771
     772      /// \brief Returns a random bool with given probability of true result.
     773      ///
     774      /// It returns a random bool with given probability of true result.
     775      bool boolean(double p) {
     776        return operator()() < p;
     777      }
     778
     779      /// Standard normal (Gauss) distribution
     780
     781      /// Standard normal (Gauss) distribution.
     782      /// \note The Cartesian form of the Box-Muller
     783      /// transformation is used to generate a random normal distribution.
     784      double gauss()
     785      {
     786        double V1,V2,S;
     787        do {
     788          V1=2*real<double>()-1;
     789          V2=2*real<double>()-1;
     790          S=V1*V1+V2*V2;
     791        } while(S>=1);
     792        return std::sqrt(-2*std::log(S)/S)*V1;
     793      }
     794      /// Normal (Gauss) distribution with given mean and standard deviation
     795
     796      /// Normal (Gauss) distribution with given mean and standard deviation.
     797      /// \sa gauss()
     798      double gauss(double mean,double std_dev)
     799      {
     800        return gauss()*std_dev+mean;
     801      }
     802
     803      /// Lognormal distribution
     804
     805      /// Lognormal distribution. The parameters are the mean and the standard
     806      /// deviation of <tt>exp(X)</tt>.
     807      ///
     808      double lognormal(double n_mean,double n_std_dev)
     809      {
     810        return std::exp(gauss(n_mean,n_std_dev));
     811      }
     812      /// Lognormal distribution
     813
     814      /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of
     815      /// the mean and the standard deviation of <tt>exp(X)</tt>.
     816      ///
     817      double lognormal(const std::pair<double,double> &params)
     818      {
     819        return std::exp(gauss(params.first,params.second));
     820      }
     821      /// Compute the lognormal parameters from mean and standard deviation
     822
     823      /// This function computes the lognormal parameters from mean and
     824      /// standard deviation. The return value can direcly be passed to
     825      /// lognormal().
     826      std::pair<double,double> lognormalParamsFromMD(double mean,
     827                                                     double std_dev)
     828      {
     829        double fr=std_dev/mean;
     830        fr*=fr;
     831        double lg=std::log(1+fr);
     832        return std::pair<double,double>(std::log(mean)-lg/2.0,std::sqrt(lg));
     833      }
     834      /// Lognormal distribution with given mean and standard deviation
     835
     836      /// Lognormal distribution with given mean and standard deviation.
     837      ///
     838      double lognormalMD(double mean,double std_dev)
     839      {
     840        return lognormal(lognormalParamsFromMD(mean,std_dev));
     841      }
     842
     843      /// Exponential distribution with given mean
     844
     845      /// This function generates an exponential distribution random number
     846      /// with mean <tt>1/lambda</tt>.
     847      ///
     848      double exponential(double lambda=1.0)
     849      {
     850        return -std::log(1.0-real<double>())/lambda;
     851      }
     852
     853      /// Gamma distribution with given integer shape
     854
     855      /// This function generates a gamma distribution random number.
     856      ///
     857      ///\param k shape parameter (<tt>k>0</tt> integer)
     858      double gamma(int k)
     859      {
     860        double s = 0;
     861        for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
     862        return s;
     863      }
     864
     865      /// Gamma distribution with given shape and scale parameter
     866
     867      /// This function generates a gamma distribution random number.
     868      ///
     869      ///\param k shape parameter (<tt>k>0</tt>)
     870      ///\param theta scale parameter
     871      ///
     872      double gamma(double k,double theta=1.0)
     873      {
     874        double xi,nu;
     875        const double delta = k-std::floor(k);
     876        const double v0=E/(E-delta);
     877        do {
     878          double V0=1.0-real<double>();
     879          double V1=1.0-real<double>();
     880          double V2=1.0-real<double>();
     881          if(V2<=v0)
     882            {
     883              xi=std::pow(V1,1.0/delta);
     884              nu=V0*std::pow(xi,delta-1.0);
     885            }
     886          else
     887            {
     888              xi=1.0-std::log(V1);
     889              nu=V0*std::exp(-xi);
     890            }
     891        } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
     892        return theta*(xi+gamma(int(std::floor(k))));
     893      }
     894
     895      /// Weibull distribution
     896
     897      /// This function generates a Weibull distribution random number.
     898      ///
     899      ///\param k shape parameter (<tt>k>0</tt>)
     900      ///\param lambda scale parameter (<tt>lambda>0</tt>)
     901      ///
     902      double weibull(double k,double lambda)
     903      {
     904        return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
     905      }
     906
     907      /// Pareto distribution
     908
     909      /// This function generates a Pareto distribution random number.
     910      ///
     911      ///\param k shape parameter (<tt>k>0</tt>)
     912      ///\param x_min location parameter (<tt>x_min>0</tt>)
     913      ///
     914      double pareto(double k,double x_min)
     915      {
     916        return exponential(gamma(k,1.0/x_min))+x_min;
     917      }
     918
     919      /// Poisson distribution
     920
     921      /// This function generates a Poisson distribution random number with
     922      /// parameter \c lambda.
     923      ///
     924      /// The probability mass function of this distribusion is
     925      /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]
     926      /// \note The algorithm is taken from the book of Donald E. Knuth titled
     927      /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the
     928      /// return value.
     929
     930      int poisson(double lambda)
     931      {
     932        const double l = std::exp(-lambda);
     933        int k=0;
     934        double p = 1.0;
     935        do {
     936          k++;
     937          p*=real<double>();
     938        } while (p>=l);
     939        return k-1;
     940      }
     941
     942      ///@}
     943
     944      ///\name Two-Dimensional Distributions
     945      ///
     946      ///@{
     947
     948      /// Uniform distribution on the full unit circle
     949
     950      /// Uniform distribution on the full unit circle.
     951      ///
     952      dim2::Point<double> disc()
     953      {
     954        double V1,V2;
     955        do {
     956          V1=2*real<double>()-1;
     957          V2=2*real<double>()-1;
     958
     959        } while(V1*V1+V2*V2>=1);
     960        return dim2::Point<double>(V1,V2);
     961      }
     962      /// A kind of two-dimensional normal (Gauss) distribution
     963
     964      /// This function provides a turning symmetric two-dimensional distribution.
     965      /// Both coordinates are of standard normal distribution, but they are not
     966      /// independent.
     967      ///
     968      /// \note The coordinates are the two random variables provided by
     969      /// the Box-Muller method.
     970      dim2::Point<double> gauss2()
     971      {
     972        double V1,V2,S;
     973        do {
     974          V1=2*real<double>()-1;
     975          V2=2*real<double>()-1;
     976          S=V1*V1+V2*V2;
     977        } while(S>=1);
     978        double W=std::sqrt(-2*std::log(S)/S);
     979        return dim2::Point<double>(W*V1,W*V2);
     980      }
     981      /// A kind of two-dimensional exponential distribution
     982
     983      /// This function provides a turning symmetric two-dimensional distribution.
     984      /// The x-coordinate is of conditionally exponential distribution
     985      /// with the condition that x is positive and y=0. If x is negative and
     986      /// y=0 then, -x is of exponential distribution. The same is true for the
     987      /// y-coordinate.
     988      dim2::Point<double> exponential2()
     989      {
     990        double V1,V2,S;
     991        do {
     992          V1=2*real<double>()-1;
     993          V2=2*real<double>()-1;
     994          S=V1*V1+V2*V2;
     995        } while(S>=1);
     996        double W=-std::log(S)/S;
     997        return dim2::Point<double>(W*V1,W*V2);
     998      }
     999
     1000      ///@}
     1001    };
     1002
     1003
     1004  };
    4691005
    4701006  /// \ingroup misc
     
    4721008  /// \brief Mersenne Twister random number generator
    4731009  ///
    474   /// The Mersenne Twister is a twisted generalized feedback
    475   /// shift-register generator of Matsumoto and Nishimura. The period
    476   /// of this generator is \f$ 2^{19937} - 1 \f$ and it is
    477   /// equi-distributed in 623 dimensions for 32-bit numbers. The time
    478   /// performance of this generator is comparable to the commonly used
    479   /// generators.
     1010  /// This class implements either the 32-bit or the 64-bit version of
     1011  /// the Mersenne Twister random number generator algorithm
     1012  /// depending on the system architecture.
     1013  ///
     1014  /// For the API description, see its base class
     1015  /// \ref _random_bits::Random.
    4801016  ///
    481   /// This implementation is specialized for both 32-bit and 64-bit
    482   /// architectures. The generators differ sligthly in the
    483   /// initialization and generation phase so they produce two
    484   /// completly different sequences.
     1017  /// \sa \ref _random_bits::Random
     1018  typedef _random_bits::Random<unsigned long> Random;
     1019
     1020  /// \ingroup misc
    4851021  ///
    486   /// The generator gives back random numbers of serveral types. To
    487   /// get a random number from a range of a floating point type you
    488   /// can use one form of the \c operator() or the \c real() member
    489   /// function. If you want to get random number from the {0, 1, ...,
    490   /// n-1} integer range use the \c operator[] or the \c integer()
    491   /// method. And to get random number from the whole range of an
    492   /// integer type you can use the argumentless \c integer() or \c
    493   /// uinteger() functions. After all you can get random bool with
    494   /// equal chance of true and false or given probability of true
    495   /// result with the \c boolean() member functions.
     1022  /// \brief Mersenne Twister random number generator (32-bit version)
    4961023  ///
    497   ///\code
    498   /// // The commented code is identical to the other
    499   /// double a = rnd();                     // [0.0, 1.0)
    500   /// // double a = rnd.real();             // [0.0, 1.0)
    501   /// double b = rnd(100.0);                // [0.0, 100.0)
    502   /// // double b = rnd.real(100.0);        // [0.0, 100.0)
    503   /// double c = rnd(1.0, 2.0);             // [1.0, 2.0)
    504   /// // double c = rnd.real(1.0, 2.0);     // [1.0, 2.0)
    505   /// int d = rnd[100000];                  // 0..99999
    506   /// // int d = rnd.integer(100000);       // 0..99999
    507   /// int e = rnd[6] + 1;                   // 1..6
    508   /// // int e = rnd.integer(1, 1 + 6);     // 1..6
    509   /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1
    510   /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1
    511   /// bool g = rnd.boolean();               // P(g = true) = 0.5
    512   /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8
    513   ///\endcode
     1024  /// This class implements the 32-bit version of the Mersenne Twister
     1025  /// random number generator algorithm. It is recommended to be used
     1026  /// when someone wants to make sure that the \e same pseudo random
     1027  /// sequence will be generated on every platfrom.
    5141028  ///
    515   /// LEMON provides a global instance of the random number
    516   /// generator which name is \ref lemon::rnd "rnd". Usually it is a
    517   /// good programming convenience to use this global generator to get
    518   /// random numbers.
    519   class Random {
    520   private:
    521 
    522     // Architecture word
    523     typedef unsigned long Word;
    524 
    525     _random_bits::RandomCore<Word> core;
    526     _random_bits::BoolProducer<Word> bool_producer;
    527 
    528 
    529   public:
    530 
    531     ///\name Initialization
    532     ///
    533     /// @{
    534 
    535     /// \brief Default constructor
    536     ///
    537     /// Constructor with constant seeding.
    538     Random() { core.initState(); }
    539 
    540     /// \brief Constructor with seed
    541     ///
    542     /// Constructor with seed. The current number type will be converted
    543     /// to the architecture word type.
    544     template <typename Number>
    545     Random(Number seed) {
    546       _random_bits::Initializer<Number, Word>::init(core, seed);
    547     }
    548 
    549     /// \brief Constructor with array seeding
    550     ///
    551     /// Constructor with array seeding. The given range should contain
    552     /// any number type and the numbers will be converted to the
    553     /// architecture word type.
    554     template <typename Iterator>
    555     Random(Iterator begin, Iterator end) {
    556       typedef typename std::iterator_traits<Iterator>::value_type Number;
    557       _random_bits::Initializer<Number, Word>::init(core, begin, end);
    558     }
    559 
    560     /// \brief Copy constructor
    561     ///
    562     /// Copy constructor. The generated sequence will be identical to
    563     /// the other sequence. It can be used to save the current state
    564     /// of the generator and later use it to generate the same
    565     /// sequence.
    566     Random(const Random& other) {
    567       core.copyState(other.core);
    568     }
    569 
    570     /// \brief Assign operator
    571     ///
    572     /// Assign operator. The generated sequence will be identical to
    573     /// the other sequence. It can be used to save the current state
    574     /// of the generator and later use it to generate the same
    575     /// sequence.
    576     Random& operator=(const Random& other) {
    577       if (&other != this) {
    578         core.copyState(other.core);
    579       }
    580       return *this;
    581     }
    582 
    583     /// \brief Seeding random sequence
    584     ///
    585     /// Seeding the random sequence. The current number type will be
    586     /// converted to the architecture word type.
    587     template <typename Number>
    588     void seed(Number seed) {
    589       _random_bits::Initializer<Number, Word>::init(core, seed);
    590     }
    591 
    592     /// \brief Seeding random sequence
    593     ///
    594     /// Seeding the random sequence. The given range should contain
    595     /// any number type and the numbers will be converted to the
    596     /// architecture word type.
    597     template <typename Iterator>
    598     void seed(Iterator begin, Iterator end) {
    599       typedef typename std::iterator_traits<Iterator>::value_type Number;
    600       _random_bits::Initializer<Number, Word>::init(core, begin, end);
    601     }
    602 
    603     /// \brief Seeding from file or from process id and time
    604     ///
    605     /// By default, this function calls the \c seedFromFile() member
    606     /// function with the <tt>/dev/urandom</tt> file. If it does not success,
    607     /// it uses the \c seedFromTime().
    608     /// \return Currently always \c true.
    609     bool seed() {
    610 #ifndef LEMON_WIN32
    611       if (seedFromFile("/dev/urandom", 0)) return true;
     1029  /// For the API description, see its base class
     1030  /// \ref _random_bits::Random.
     1031  ///
     1032  /// \sa \ref _random_bits::Random
     1033  typedef _random_bits::Random<unsigned int> Random32;
     1034
     1035  /// \ingroup misc
     1036  ///
     1037  /// \brief Mersenne Twister random number generator (64-bit version)
     1038  ///
     1039  /// This class implements the 64-bit version of the Mersenne Twister
     1040  /// random number generator algorithm. (Even though it runs
     1041  /// on 32-bit architectures, too.) It is recommended to be used when
     1042  /// someone wants to make sure that the \e same pseudo random sequence
     1043  /// will be generated on every platfrom.
     1044  ///
     1045  /// For the API description, see its base class
     1046  /// \ref _random_bits::Random.
     1047  ///
     1048  /// \sa \ref _random_bits::Random
     1049  typedef _random_bits::Random<unsigned long long> Random64;
     1050
     1051  extern Random rnd;
     1052 
     1053}
     1054
    6121055#endif
    613       if (seedFromTime()) return true;
    614       return false;
    615     }
    616 
    617     /// \brief Seeding from file
    618     ///
    619     /// Seeding the random sequence from file. The linux kernel has two
    620     /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which
    621     /// could give good seed values for pseudo random generators (The
    622     /// difference between two devices is that the <tt>random</tt> may
    623     /// block the reading operation while the kernel can give good
    624     /// source of randomness, while the <tt>urandom</tt> does not
    625     /// block the input, but it could give back bytes with worse
    626     /// entropy).
    627     /// \param file The source file
    628     /// \param offset The offset, from the file read.
    629     /// \return \c true when the seeding successes.
    630 #ifndef LEMON_WIN32
    631     bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)
    632 #else
    633     bool seedFromFile(const std::string& file = "", int offset = 0)
    634 #endif
    635     {
    636       std::ifstream rs(file.c_str());
    637       const int size = 4;
    638       Word buf[size];
    639       if (offset != 0 && !rs.seekg(offset)) return false;
    640       if (!rs.read(reinterpret_cast<char*>(buf), sizeof(buf))) return false;
    641       seed(buf, buf + size);
    642       return true;
    643     }
    644 
    645     /// \brief Seding from process id and time
    646     ///
    647     /// Seding from process id and time. This function uses the
    648     /// current process id and the current time for initialize the
    649     /// random sequence.
    650     /// \return Currently always \c true.
    651     bool seedFromTime() {
    652 #ifndef LEMON_WIN32
    653       timeval tv;
    654       gettimeofday(&tv, 0);
    655       seed(getpid() + tv.tv_sec + tv.tv_usec);
    656 #else
    657       seed(bits::getWinRndSeed());
    658 #endif
    659       return true;
    660     }
    661 
    662     /// @}
    663 
    664     ///\name Uniform Distributions
    665     ///
    666     /// @{
    667 
    668     /// \brief Returns a random real number from the range [0, 1)
    669     ///
    670     /// It returns a random real number from the range [0, 1). The
    671     /// default Number type is \c double.
    672     template <typename Number>
    673     Number real() {
    674       return _random_bits::RealConversion<Number, Word>::convert(core);
    675     }
    676 
    677     double real() {
    678       return real<double>();
    679     }
    680 
    681     /// \brief Returns a random real number from the range [0, 1)
    682     ///
    683     /// It returns a random double from the range [0, 1).
    684     double operator()() {
    685       return real<double>();
    686     }
    687 
    688     /// \brief Returns a random real number from the range [0, b)
    689     ///
    690     /// It returns a random real number from the range [0, b).
    691     double operator()(double b) {
    692       return real<double>() * b;
    693     }
    694 
    695     /// \brief Returns a random real number from the range [a, b)
    696     ///
    697     /// It returns a random real number from the range [a, b).
    698     double operator()(double a, double b) {
    699       return real<double>() * (b - a) + a;
    700     }
    701 
    702     /// \brief Returns a random integer from a range
    703     ///
    704     /// It returns a random integer from the range {0, 1, ..., b - 1}.
    705     template <typename Number>
    706     Number integer(Number b) {
    707       return _random_bits::Mapping<Number, Word>::map(core, b);
    708     }
    709 
    710     /// \brief Returns a random integer from a range
    711     ///
    712     /// It returns a random integer from the range {a, a + 1, ..., b - 1}.
    713     template <typename Number>
    714     Number integer(Number a, Number b) {
    715       return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;
    716     }
    717 
    718     /// \brief Returns a random integer from a range
    719     ///
    720     /// It returns a random integer from the range {0, 1, ..., b - 1}.
    721     template <typename Number>
    722     Number operator[](Number b) {
    723       return _random_bits::Mapping<Number, Word>::map(core, b);
    724     }
    725 
    726     /// \brief Returns a random non-negative integer
    727     ///
    728     /// It returns a random non-negative integer uniformly from the
    729     /// whole range of the current \c Number type. The default result
    730     /// type of this function is <tt>unsigned int</tt>.
    731     template <typename Number>
    732     Number uinteger() {
    733       return _random_bits::IntConversion<Number, Word>::convert(core);
    734     }
    735 
    736     unsigned int uinteger() {
    737       return uinteger<unsigned int>();
    738     }
    739 
    740     /// \brief Returns a random integer
    741     ///
    742     /// It returns a random integer uniformly from the whole range of
    743     /// the current \c Number type. The default result type of this
    744     /// function is \c int.
    745     template <typename Number>
    746     Number integer() {
    747       static const int nb = std::numeric_limits<Number>::digits +
    748         (std::numeric_limits<Number>::is_signed ? 1 : 0);
    749       return _random_bits::IntConversion<Number, Word, nb>::convert(core);
    750     }
    751 
    752     int integer() {
    753       return integer<int>();
    754     }
    755 
    756     /// \brief Returns a random bool
    757     ///
    758     /// It returns a random bool. The generator holds a buffer for
    759     /// random bits. Every time when it become empty the generator makes
    760     /// a new random word and fill the buffer up.
    761     bool boolean() {
    762       return bool_producer.convert(core);
    763     }
    764 
    765     /// @}
    766 
    767     ///\name Non-uniform Distributions
    768     ///
    769     ///@{
    770 
    771     /// \brief Returns a random bool with given probability of true result.
    772     ///
    773     /// It returns a random bool with given probability of true result.
    774     bool boolean(double p) {
    775       return operator()() < p;
    776     }
    777 
    778     /// Standard normal (Gauss) distribution
    779 
    780     /// Standard normal (Gauss) distribution.
    781     /// \note The Cartesian form of the Box-Muller
    782     /// transformation is used to generate a random normal distribution.
    783     double gauss()
    784     {
    785       double V1,V2,S;
    786       do {
    787         V1=2*real<double>()-1;
    788         V2=2*real<double>()-1;
    789         S=V1*V1+V2*V2;
    790       } while(S>=1);
    791       return std::sqrt(-2*std::log(S)/S)*V1;
    792     }
    793     /// Normal (Gauss) distribution with given mean and standard deviation
    794 
    795     /// Normal (Gauss) distribution with given mean and standard deviation.
    796     /// \sa gauss()
    797     double gauss(double mean,double std_dev)
    798     {
    799       return gauss()*std_dev+mean;
    800     }
    801 
    802     /// Lognormal distribution
    803 
    804     /// Lognormal distribution. The parameters are the mean and the standard
    805     /// deviation of <tt>exp(X)</tt>.
    806     ///
    807     double lognormal(double n_mean,double n_std_dev)
    808     {
    809       return std::exp(gauss(n_mean,n_std_dev));
    810     }
    811     /// Lognormal distribution
    812 
    813     /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of
    814     /// the mean and the standard deviation of <tt>exp(X)</tt>.
    815     ///
    816     double lognormal(const std::pair<double,double> &params)
    817     {
    818       return std::exp(gauss(params.first,params.second));
    819     }
    820     /// Compute the lognormal parameters from mean and standard deviation
    821 
    822     /// This function computes the lognormal parameters from mean and
    823     /// standard deviation. The return value can direcly be passed to
    824     /// lognormal().
    825     std::pair<double,double> lognormalParamsFromMD(double mean,
    826                                                    double std_dev)
    827     {
    828       double fr=std_dev/mean;
    829       fr*=fr;
    830       double lg=std::log(1+fr);
    831       return std::pair<double,double>(std::log(mean)-lg/2.0,std::sqrt(lg));
    832     }
    833     /// Lognormal distribution with given mean and standard deviation
    834 
    835     /// Lognormal distribution with given mean and standard deviation.
    836     ///
    837     double lognormalMD(double mean,double std_dev)
    838     {
    839       return lognormal(lognormalParamsFromMD(mean,std_dev));
    840     }
    841 
    842     /// Exponential distribution with given mean
    843 
    844     /// This function generates an exponential distribution random number
    845     /// with mean <tt>1/lambda</tt>.
    846     ///
    847     double exponential(double lambda=1.0)
    848     {
    849       return -std::log(1.0-real<double>())/lambda;
    850     }
    851 
    852     /// Gamma distribution with given integer shape
    853 
    854     /// This function generates a gamma distribution random number.
    855     ///
    856     ///\param k shape parameter (<tt>k>0</tt> integer)
    857     double gamma(int k)
    858     {
    859       double s = 0;
    860       for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
    861       return s;
    862     }
    863 
    864     /// Gamma distribution with given shape and scale parameter
    865 
    866     /// This function generates a gamma distribution random number.
    867     ///
    868     ///\param k shape parameter (<tt>k>0</tt>)
    869     ///\param theta scale parameter
    870     ///
    871     double gamma(double k,double theta=1.0)
    872     {
    873       double xi,nu;
    874       const double delta = k-std::floor(k);
    875       const double v0=E/(E-delta);
    876       do {
    877         double V0=1.0-real<double>();
    878         double V1=1.0-real<double>();
    879         double V2=1.0-real<double>();
    880         if(V2<=v0)
    881           {
    882             xi=std::pow(V1,1.0/delta);
    883             nu=V0*std::pow(xi,delta-1.0);
    884           }
    885         else
    886           {
    887             xi=1.0-std::log(V1);
    888             nu=V0*std::exp(-xi);
    889           }
    890       } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
    891       return theta*(xi+gamma(int(std::floor(k))));
    892     }
    893 
    894     /// Weibull distribution
    895 
    896     /// This function generates a Weibull distribution random number.
    897     ///
    898     ///\param k shape parameter (<tt>k>0</tt>)
    899     ///\param lambda scale parameter (<tt>lambda>0</tt>)
    900     ///
    901     double weibull(double k,double lambda)
    902     {
    903       return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
    904     }
    905 
    906     /// Pareto distribution
    907 
    908     /// This function generates a Pareto distribution random number.
    909     ///
    910     ///\param k shape parameter (<tt>k>0</tt>)
    911     ///\param x_min location parameter (<tt>x_min>0</tt>)
    912     ///
    913     double pareto(double k,double x_min)
    914     {
    915       return exponential(gamma(k,1.0/x_min))+x_min;
    916     }
    917 
    918     /// Poisson distribution
    919 
    920     /// This function generates a Poisson distribution random number with
    921     /// parameter \c lambda.
    922     ///
    923     /// The probability mass function of this distribusion is
    924     /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]
    925     /// \note The algorithm is taken from the book of Donald E. Knuth titled
    926     /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the
    927     /// return value.
    928 
    929     int poisson(double lambda)
    930     {
    931       const double l = std::exp(-lambda);
    932       int k=0;
    933       double p = 1.0;
    934       do {
    935         k++;
    936         p*=real<double>();
    937       } while (p>=l);
    938       return k-1;
    939     }
    940 
    941     ///@}
    942 
    943     ///\name Two Dimensional Distributions
    944     ///
    945     ///@{
    946 
    947     /// Uniform distribution on the full unit circle
    948 
    949     /// Uniform distribution on the full unit circle.
    950     ///
    951     dim2::Point<double> disc()
    952     {
    953       double V1,V2;
    954       do {
    955         V1=2*real<double>()-1;
    956         V2=2*real<double>()-1;
    957 
    958       } while(V1*V1+V2*V2>=1);
    959       return dim2::Point<double>(V1,V2);
    960     }
    961     /// A kind of two dimensional normal (Gauss) distribution
    962 
    963     /// This function provides a turning symmetric two-dimensional distribution.
    964     /// Both coordinates are of standard normal distribution, but they are not
    965     /// independent.
    966     ///
    967     /// \note The coordinates are the two random variables provided by
    968     /// the Box-Muller method.
    969     dim2::Point<double> gauss2()
    970     {
    971       double V1,V2,S;
    972       do {
    973         V1=2*real<double>()-1;
    974         V2=2*real<double>()-1;
    975         S=V1*V1+V2*V2;
    976       } while(S>=1);
    977       double W=std::sqrt(-2*std::log(S)/S);
    978       return dim2::Point<double>(W*V1,W*V2);
    979     }
    980     /// A kind of two dimensional exponential distribution
    981 
    982     /// This function provides a turning symmetric two-dimensional distribution.
    983     /// The x-coordinate is of conditionally exponential distribution
    984     /// with the condition that x is positive and y=0. If x is negative and
    985     /// y=0 then, -x is of exponential distribution. The same is true for the
    986     /// y-coordinate.
    987     dim2::Point<double> exponential2()
    988     {
    989       double V1,V2,S;
    990       do {
    991         V1=2*real<double>()-1;
    992         V2=2*real<double>()-1;
    993         S=V1*V1+V2*V2;
    994       } while(S>=1);
    995       double W=-std::log(S)/S;
    996       return dim2::Point<double>(W*V1,W*V2);
    997     }
    998 
    999     ///@}
    1000   };
    1001 
    1002 
    1003   extern Random rnd;
    1004 
    1005 }
    1006 
    1007 #endif
Note: See TracChangeset for help on using the changeset viewer.