463 --num; |
463 --num; |
464 return r; |
464 return r; |
465 } |
465 } |
466 }; |
466 }; |
467 |
467 |
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> ¶ms) |
|
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 }; |
469 |
1006 |
470 /// \ingroup misc |
1007 /// \ingroup misc |
471 /// |
1008 /// |
472 /// \brief Mersenne Twister random number generator |
1009 /// \brief Mersenne Twister random number generator |
473 /// |
1010 /// |
474 /// The Mersenne Twister is a twisted generalized feedback |
1011 /// This class implements either the 32 bit or the 64 bit version of |
475 /// shift-register generator of Matsumoto and Nishimura. The period |
1012 /// the Mersenne Twister random number generator algorithm |
476 /// of this generator is \f$ 2^{19937} - 1 \f$ and it is |
1013 /// depending the the system architecture. |
477 /// equi-distributed in 623 dimensions for 32-bit numbers. The time |
1014 /// |
478 /// performance of this generator is comparable to the commonly used |
1015 /// For the API description, see its base class \ref |
479 /// generators. |
1016 /// _random_bits::Random |
480 /// |
1017 /// |
481 /// This implementation is specialized for both 32-bit and 64-bit |
1018 /// \sa \ref _random_bits::Random |
482 /// architectures. The generators differ sligthly in the |
1019 typedef _random_bits::Random<unsigned long> Random; |
483 /// initialization and generation phase so they produce two |
1020 /// \ingroup misc |
484 /// completly different sequences. |
|
485 /// |
1021 /// |
486 /// The generator gives back random numbers of serveral types. To |
1022 /// \brief Mersenne Twister random number generator (32 bit version) |
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. |
|
496 /// |
1023 /// |
497 ///\code |
1024 /// This class implements the 32 bit version of the Mersenne Twister |
498 /// // The commented code is identical to the other |
1025 /// random number generator algorithm. It is recommended to be used |
499 /// double a = rnd(); // [0.0, 1.0) |
1026 /// when someone wants to make sure that the \e same pseudo random |
500 /// // double a = rnd.real(); // [0.0, 1.0) |
1027 /// sequence will be generated on every platfrom. |
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 |
|
514 /// |
1028 /// |
515 /// LEMON provides a global instance of the random number |
1029 /// For the API description, see its base class \ref |
516 /// generator which name is \ref lemon::rnd "rnd". Usually it is a |
1030 /// _random_bits::Random |
517 /// good programming convenience to use this global generator to get |
1031 /// |
518 /// random numbers. |
1032 /// \sa \ref _random_bits::Random |
519 class Random { |
1033 |
520 private: |
1034 typedef _random_bits::Random<unsigned int> Random32; |
521 |
1035 /// \ingroup misc |
522 // Architecture word |
1036 /// |
523 typedef unsigned long Word; |
1037 /// \brief Mersenne Twister random number generator (64 bit version) |
524 |
1038 /// |
525 _random_bits::RandomCore<Word> core; |
1039 /// This class implements the 64 bit version of the Mersenne Twister |
526 _random_bits::BoolProducer<Word> bool_producer; |
1040 /// random number generator algorithm. (Even though it will run |
527 |
1041 /// on 32 bit architectures, too.) It is recommended to ber used when |
528 |
1042 /// someone wants to make sure that the \e same pseudo random sequence |
529 public: |
1043 /// will be generated on every platfrom. |
530 |
1044 /// |
531 ///\name Initialization |
1045 /// For the API description, see its base class \ref |
532 /// |
1046 /// _random_bits::Random |
533 /// @{ |
1047 /// |
534 |
1048 /// \sa \ref _random_bits::Random |
535 /// \brief Default constructor |
1049 typedef _random_bits::Random<unsigned long long> Random64; |
536 /// |
1050 |
537 /// Constructor with constant seeding. |
1051 |
538 Random() { core.initState(); } |
1052 extern Random rnd; |
539 |
1053 |
540 /// \brief Constructor with seed |
1054 |
541 /// |
1055 } |
542 /// Constructor with seed. The current number type will be converted |
1056 |
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; |
|
612 #endif |
1057 #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> ¶ms) |
|
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 |
|