Changeset 1404:c8d0179a32a2 in lemon for lemon/random.h
- Timestamp:
- 10/17/18 19:22:52 (6 years ago)
- 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
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
lemon/random.h
r1380 r1404 340 340 num = rnd() & mask; 341 341 } while (num > max); 342 return num;342 return static_cast<Result>(num); 343 343 } 344 344 }; -
lemon/random.h
r1396 r1404 111 111 static const Word loMask = (1u << 31) - 1; 112 112 static const Word hiMask = ~loMask; 113 114 113 115 114 static Word tempering(Word rnd) { … … 244 243 private: 245 244 246 247 245 void fillState() { 248 246 static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask }; … … 252 250 current = state + length; 253 251 254 registerWord *curr = state + length - 1;255 registerlong num;252 Word *curr = state + length - 1; 253 long num; 256 254 257 255 num = length - shift; … … 272 270 } 273 271 274 275 272 Word *current; 276 273 Word state[length]; … … 297 294 template <typename Result, typename Word, 298 295 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)> 300 297 struct IntConversion { 301 298 static const int bits = std::numeric_limits<Word>::digits; … … 466 463 }; 467 464 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> ¶ms) 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 }; 469 1005 470 1006 /// \ingroup misc … … 472 1008 /// \brief Mersenne Twister random number generator 473 1009 /// 474 /// Th e Mersenne Twister is a twisted generalized feedback475 /// shift-register generator of Matsumoto and Nishimura. The period476 /// of this generator is \f$ 2^{19937} - 1 \f$ and it is477 /// equi-distributed in 623 dimensions for 32-bit numbers. The time478 /// performance of this generator is comparable to the commonly used479 /// 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. 480 1016 /// 481 /// This implementation is specialized for both 32-bit and 64-bit482 /// architectures. The generators differ sligthly in the483 /// 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 485 1021 /// 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) 496 1023 /// 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. 514 1028 /// 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 612 1055 #endif 613 if (seedFromTime()) return true;614 return false;615 }616 617 /// \brief Seeding from file618 ///619 /// Seeding the random sequence from file. The linux kernel has two620 /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which621 /// could give good seed values for pseudo random generators (The622 /// difference between two devices is that the <tt>random</tt> may623 /// block the reading operation while the kernel can give good624 /// source of randomness, while the <tt>urandom</tt> does not625 /// block the input, but it could give back bytes with worse626 /// entropy).627 /// \param file The source file628 /// \param offset The offset, from the file read.629 /// \return \c true when the seeding successes.630 #ifndef LEMON_WIN32631 bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)632 #else633 bool seedFromFile(const std::string& file = "", int offset = 0)634 #endif635 {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 time646 ///647 /// Seding from process id and time. This function uses the648 /// current process id and the current time for initialize the649 /// random sequence.650 /// \return Currently always \c true.651 bool seedFromTime() {652 #ifndef LEMON_WIN32653 timeval tv;654 gettimeofday(&tv, 0);655 seed(getpid() + tv.tv_sec + tv.tv_usec);656 #else657 seed(bits::getWinRndSeed());658 #endif659 return true;660 }661 662 /// @}663 664 ///\name Uniform Distributions665 ///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). The671 /// 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 range703 ///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 range711 ///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 range719 ///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 integer727 ///728 /// It returns a random non-negative integer uniformly from the729 /// whole range of the current \c Number type. The default result730 /// 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 integer741 ///742 /// It returns a random integer uniformly from the whole range of743 /// the current \c Number type. The default result type of this744 /// 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 bool757 ///758 /// It returns a random bool. The generator holds a buffer for759 /// random bits. Every time when it become empty the generator makes760 /// 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 Distributions768 ///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) distribution779 780 /// Standard normal (Gauss) distribution.781 /// \note The Cartesian form of the Box-Muller782 /// 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 deviation794 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 distribution803 804 /// Lognormal distribution. The parameters are the mean and the standard805 /// 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 distribution812 813 /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of814 /// the mean and the standard deviation of <tt>exp(X)</tt>.815 ///816 double lognormal(const std::pair<double,double> ¶ms)817 {818 return std::exp(gauss(params.first,params.second));819 }820 /// Compute the lognormal parameters from mean and standard deviation821 822 /// This function computes the lognormal parameters from mean and823 /// standard deviation. The return value can direcly be passed to824 /// 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 deviation834 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 mean843 844 /// This function generates an exponential distribution random number845 /// 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 shape853 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 parameter865 866 /// This function generates a gamma distribution random number.867 ///868 ///\param k shape parameter (<tt>k>0</tt>)869 ///\param theta scale parameter870 ///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 else886 {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 distribution895 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 distribution907 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 distribution919 920 /// This function generates a Poisson distribution random number with921 /// parameter \c lambda.922 ///923 /// The probability mass function of this distribusion is924 /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]925 /// \note The algorithm is taken from the book of Donald E. Knuth titled926 /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the927 /// 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 Distributions944 ///945 ///@{946 947 /// Uniform distribution on the full unit circle948 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) distribution962 963 /// This function provides a turning symmetric two-dimensional distribution.964 /// Both coordinates are of standard normal distribution, but they are not965 /// independent.966 ///967 /// \note The coordinates are the two random variables provided by968 /// 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 distribution981 982 /// This function provides a turning symmetric two-dimensional distribution.983 /// The x-coordinate is of conditionally exponential distribution984 /// with the condition that x is positive and y=0. If x is negative and985 /// y=0 then, -x is of exponential distribution. The same is true for the986 /// 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.