COIN-OR::LEMON - Graph Library

Changeset 1379:db1d342a1087 in lemon


Ignore:
Timestamp:
10/08/15 13:48:09 (9 years ago)
Author:
Alpar Juttner <alpar@…>
Branch:
default
Phase:
public
Message:

Platform independent Random generators (#602)

Files:
2 edited

Legend:

Unmodified
Added
Removed
  • lemon/random.h

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

    r463 r1379  
    2222int seed_array[] = {1, 2};
    2323
     24int rnd_seq32[] = {
     252732, 43567, 42613, 52416, 45891, 21243, 30403, 32103,
     2662501, 33003, 12172, 5192, 32511, 50057, 43723, 7813,
     2723720, 35343, 6637, 30280, 44566, 31019, 18898, 33867,
     285994, 1688, 11513, 59011, 48056, 25544, 39168, 25365,
     2917530, 8366, 27063, 49861, 55169, 63848, 11863, 49608
     30};
     31int rnd_seq64[] = {
     3256382, 63883, 59577, 64750, 9644, 59886, 57647, 18152,
     3328520, 64078, 17818, 49294, 26424, 26697, 53684, 19209,
     3435404, 12121, 12837, 11827, 32156, 58333, 62553, 7907,
     3564427, 39399, 21971, 48789, 46981, 15716, 53335, 65256,
     3612999, 15308, 10906, 42162, 47587, 43006, 53921, 18716
     37};
     38
     39void seq_test() {
     40  for(int i=0;i<5;i++) {
     41    lemon::Random32 r32(i);
     42    lemon::Random64 r64(i);
     43    for(int j=0;j<8;j++) {
     44      check(r32[65536]==rnd_seq32[i*8+j], "Wrong random sequence");
     45      check(r64[65536]==rnd_seq64[i*8+j], "Wrong random sequence");
     46    }
     47  }
     48}
     49
     50
    2451int main()
    2552{
     
    3764                  (sizeof(seed_array) / sizeof(seed_array[0])));
    3865
     66  seq_test();
    3967  return 0;
    4068}
Note: See TracChangeset for help on using the changeset viewer.