RhodeCode

/* -*- C++ -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library

 *

 * Copyright (C) 2003-2007

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

///\file

///\brief Instantiation of the Random class.

#include <lemon/random.h>

namespace lemon {

  /// \brief Global random number generator instance

  ///

  Random rnd;

}

 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT

 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS

 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE

 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,

 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES

 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR

 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)

 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,

 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)

 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED

 * OF THE POSSIBILITY OF SUCH DAMAGE.

 *

 *

 * Any feedback is very welcome.

 * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html

 * email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)

 */

#ifndef LEMON_RANDOM_H

#define LEMON_RANDOM_H

#include <algorithm>

#include <iterator>

#include <vector>

#include <ctime>

#include <cmath>

#include <lemon/dim2.h>

///\ingroup misc

///\file

///\brief Mersenne Twister random number generator

namespace lemon {

  namespace _random_bits {

    template <typename _Word, int _bits = std::numeric_limits<_Word>::digits>

    struct RandomTraits {};

    template <typename _Word>

    struct RandomTraits<_Word, 32> {

      typedef _Word Word;

      static const int bits = 32;

      static const int length = 624;

      static const int shift = 397;

      static const Word mul = 0x6c078965u;

      static const Word arrayInit = 0x012BD6AAu;

      static const Word arrayMul1 = 0x0019660Du;

      static const Word arrayMul2 = 0x5D588B65u;

      static const Word mask = 0x9908B0DFu;

      static const Word loMask = (1u << 31) - 1;

      static const Word hiMask = ~loMask;

      static Word tempering(Word rnd) {

        rnd ^= (rnd >> 11);

        rnd ^= (rnd << 7) & 0x9D2C5680u;

        rnd ^= (rnd << 15) & 0xEFC60000u;

        rnd ^= (rnd >> 18);

  /// get a random number from a range of a floating point type you

  /// can use one form of the \c operator() or the \c real() member

  /// function. If you want to get random number from the {0, 1, ...,

  /// n-1} integer range use the \c operator[] or the \c integer()

  /// method. And to get random number from the whole range of an

  /// integer type you can use the argumentless \c integer() or \c

  /// uinteger() functions. After all you can get random bool with

  /// equal chance of true and false or given probability of true

  /// result with the \c boolean() member functions.

  ///

  ///\code

  /// // The commented code is identical to the other

  /// double a = rnd();                     // [0.0, 1.0)

  /// // double a = rnd.real();             // [0.0, 1.0)

  /// double b = rnd(100.0);                // [0.0, 100.0)

  /// // double b = rnd.real(100.0);        // [0.0, 100.0)

  /// double c = rnd(1.0, 2.0);             // [1.0, 2.0)

  /// // double c = rnd.real(1.0, 2.0);     // [1.0, 2.0)

  /// int d = rnd[100000];                  // 0..99999

  /// // int d = rnd.integer(100000);       // 0..99999

  /// int e = rnd[6] + 1;                   // 1..6

  /// // int e = rnd.integer(1, 1 + 6);     // 1..6

  /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1

  /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1

  /// bool g = rnd.boolean();               // P(g = true) = 0.5

  /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8

  ///\endcode

  ///

  /// The lemon provides a global instance of the random number

  /// generator which name is \ref lemon::rnd "rnd". Usually it is a

  /// good programming convenience to use this global generator to get

  /// random numbers.

  class Random {

  private:

    typedef unsigned long Word;

    _random_bits::RandomCore<Word> core;

    _random_bits::BoolProducer<Word> bool_producer;

  public:

    /// \brief Constructor

    ///

    /// Constructor with constant seeding.

    Random() { core.initState(); }

    /// \brief Constructor

    ///

    /// Constructor with seed. The current number type will be converted

    /// to the architecture word type.

    template <typename Number>

    Random(Number seed) { 

      _random_bits::Initializer<Number, Word>::init(core, seed);

    }

    /// \brief Constructor

    ///

    /// Constructor with array seeding. The given range should contain

    /// any number type and the numbers will be converted to the

    /// architecture word type.

    template <typename Iterator>

    Random(Iterator begin, Iterator end) { 

      typedef typename std::iterator_traits<Iterator>::value_type Number;

      _random_bits::Initializer<Number, Word>::init(core, begin, end);

    }

    }

    ///\name Nonuniform distributions

    ///

    ///@{

    /// \brief Returns a random bool

    ///

    /// It returns a random bool with given probability of true result

    bool boolean(double p) {

      return operator()() < p;

    }

    /// Standard Gauss distribution

    /// Standard Gauss distribution.

    /// \note The Cartesian form of the Box-Muller

    /// transformation is used to generate a random normal distribution.

    /// \todo Consider using the "ziggurat" method instead.

    double gauss() 

    {

      double V1,V2,S;

      do {

	V1=2*real<double>()-1;

	V2=2*real<double>()-1;

	S=V1*V1+V2*V2;

      } while(S>=1);

      return std::sqrt(-2*std::log(S)/S)*V1;

    }

    /// Gauss distribution with given mean and standard deviation

    /// \sa gauss()

    double gauss(double mean,double std_dev)

    {

      return gauss()*std_dev+mean;

    }

    /// Exponential distribution with given mean

    /// This function generates an exponential distribution random number

    /// with mean <tt>1/lambda</tt>.

    ///

    double exponential(double lambda=1.0)

    {

      return -std::log(1.0-real<double>())/lambda;

    }

    /// Gamma distribution with given integer shape

    /// This function generates a gamma distribution random number.

    /// 

    ///\param k shape parameter (<tt>k>0</tt> integer)

    double gamma(int k) 

    {

      double s = 0;

      for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());

      return s;

    }

    /// Gamma distribution with given shape and scale parameter

    /// This function generates a gamma distribution random number.

    /// 

    ///\param k shape parameter (<tt>k>0</tt>)

    /// Weibull distribution

    /// This function generates a Weibull distribution random number.

    /// 

    ///\param k shape parameter (<tt>k>0</tt>)

    ///\param lambda scale parameter (<tt>lambda>0</tt>)

    ///

    double weibull(double k,double lambda)

    {

      return lambda*pow(-std::log(1.0-real<double>()),1.0/k);

    }  

    /// Pareto distribution

    /// This function generates a Pareto distribution random number.

    /// 

    ///\param k shape parameter (<tt>k>0</tt>)

    ///\param x_min location parameter (<tt>x_min>0</tt>)

    ///

    double pareto(double k,double x_min)

    {

      return exponential(gamma(k,1.0/x_min));

    }  

    ///@}

    ///\name Two dimensional distributions

    ///

    ///@{

    /// Uniform distribution on the full unit circle.

    dim2::Point<double> disc() 

    {

      double V1,V2;

      do {

	V1=2*real<double>()-1;

	V2=2*real<double>()-1;

      } while(V1*V1+V2*V2>=1);

      return dim2::Point<double>(V1,V2);

    }

    /// A kind of two dimensional Gauss distribution

    /// This function provides a turning symmetric two-dimensional distribution.

    /// Both coordinates are of standard normal distribution, but they are not

    /// independent.

    ///

    /// \note The coordinates are the two random variables provided by

    /// the Box-Muller method.

    dim2::Point<double> gauss2()

    {

      double V1,V2,S;

      do {

	V1=2*real<double>()-1;

	V2=2*real<double>()-1;

	S=V1*V1+V2*V2;

      } while(S>=1);

      double W=std::sqrt(-2*std::log(S)/S);

      return dim2::Point<double>(W*V1,W*V2);

    }

    /// A kind of two dimensional exponential distribution

    /// This function provides a turning symmetric two-dimensional distribution.

#ifndef LEMON_TOLERANCE_H

#define LEMON_TOLERANCE_H

///\ingroup misc

///\file

///\brief A basic tool to handle the anomalies of calculation with

///floating point numbers.

///

///\todo It should be in a module like "Basic tools"

namespace lemon {

  /// \addtogroup misc

  /// @{

  ///\brief A class to provide a basic way to

  ///handle the comparison of numbers that are obtained

  ///as a result of a probably inexact computation.

  ///

  ///Tolerance is a class to provide a basic way to

  ///handle the comparison of numbers that are obtained

  ///as a result of a probably inexact computation.

  ///

  ///This is an abstract class, it should be specialized for all numerical

  ///data types. These specialized classes like \ref Tolerance\<double\>

  ///may offer additional tuning parameters.

  ///

  ///\sa Tolerance<float>

  ///\sa Tolerance<double>

  ///\sa Tolerance<long double>

  ///\sa Tolerance<int>

  ///\sa Tolerance<long long int>

  ///\sa Tolerance<unsigned int>

  ///\sa Tolerance<unsigned long long int>

  template<class T>

  class Tolerance

  {

  public:

    typedef T Value;

    ///\name Comparisons

    ///The concept is that these bool functions return with \c true only if

    ///the related comparisons hold even if some numerical error appeared

    ///during the computations.

    ///@{

    ///Returns \c true if \c a is \e surely strictly less than \c b

    static bool less(Value a,Value b) {return false;}

    ///Returns \c true if \c a is \e surely different from \c b

    static bool different(Value a,Value b) {return false;}

    ///Returns \c true if \c a is \e surely positive

    static bool positive(Value a) {return false;}

    ///Returns \c true if \c a is \e surely negative

    static bool negative(Value a) {return false;}

    ///Returns \c true if \c a is \e surely non-zero

    static bool nonZero(Value a) {return false;}

    ///@}

    ///Returns the zero value.

    static Value zero() {return T();}

    //   static bool finite(Value a) {}

    //   static Value big() {}

    ///\e

    typedef float Value;

    ///Constructor setting the epsilon tolerance to the default value.

    Tolerance() : _epsilon(def_epsilon) {}

    ///Constructor setting the epsilon tolerance.

    Tolerance(float e) : _epsilon(e) {}

    ///Return the epsilon value.

    Value epsilon() const {return _epsilon;}

    ///Set the epsilon value.

    void epsilon(Value e) {_epsilon=e;}

    ///Return the default epsilon value.

    static Value defaultEpsilon() {return def_epsilon;}

    ///Set the default epsilon value.

    static void defaultEpsilon(Value e) {def_epsilon=e;}

    ///\name Comparisons

    ///See class Tolerance for more details.

    ///@{

    ///Returns \c true if \c a is \e surely strictly less than \c b

    bool less(Value a,Value b) const {return a+_epsilon<b;}

    ///Returns \c true if \c a is \e surely different from \c b

    bool different(Value a,Value b) const { return less(a,b)||less(b,a); }

    ///Returns \c true if \c a is \e surely positive

    bool positive(Value a) const { return _epsilon<a; }

    ///Returns \c true if \c a is \e surely negative

    bool negative(Value a) const { return -_epsilon>a; }

    ///Returns \c true if \c a is \e surely non-zero

    ///@}

    ///Returns zero

    static Value zero() {return 0;}

  };

  ///Double specialization of \ref Tolerance.

  ///Double specialization of \ref Tolerance.

  ///\sa Tolerance

  ///\relates Tolerance

  template<>

  class Tolerance<double>

  {

    static double def_epsilon;

    double _epsilon;

  public:

    ///\e

    typedef double Value;

    ///Constructor setting the epsilon tolerance to the default value.

    Tolerance() : _epsilon(def_epsilon) {}

    ///Constructor setting the epsilon tolerance.

    Tolerance(double e) : _epsilon(e) {}

    ///Return the epsilon value.

    Value epsilon() const {return _epsilon;}

    ///Set the epsilon value.

    void epsilon(Value e) {_epsilon=e;}

    ///Return the default epsilon value.

    static Value defaultEpsilon() {return def_epsilon;}

    ///Set the default epsilon value.

    static void defaultEpsilon(Value e) {def_epsilon=e;}

    ///\name Comparisons

    ///See class Tolerance for more details.

    ///@{

    ///Returns \c true if \c a is \e surely strictly less than \c b

    bool less(Value a,Value b) const {return a+_epsilon<b;}

    ///Returns \c true if \c a is \e surely different from \c b

    bool different(Value a,Value b) const { return less(a,b)||less(b,a); }

    ///Returns \c true if \c a is \e surely positive

    bool positive(Value a) const { return _epsilon<a; }

    ///Returns \c true if \c a is \e surely negative

    bool negative(Value a) const { return -_epsilon>a; }

    ///Returns \c true if \c a is \e surely non-zero

    ///@}

    ///Returns zero

    static Value zero() {return 0;}

  };

  ///Long double specialization of \ref Tolerance.

  ///Long double specialization of \ref Tolerance.

  ///\sa Tolerance

  ///\relates Tolerance

  template<>

  class Tolerance<long double>

  {

    static long double def_epsilon;

    long double _epsilon;

  public:

    ///\e

    typedef long double Value;

    ///Constructor setting the epsilon tolerance to the default value.

    Tolerance() : _epsilon(def_epsilon) {}

    ///Constructor setting the epsilon tolerance.

    Tolerance(long double e) : _epsilon(e) {}

    ///Return the epsilon value.

    Value epsilon() const {return _epsilon;}

    ///Set the epsilon value.

    void epsilon(Value e) {_epsilon=e;}

    ///Return the default epsilon value.

    static Value defaultEpsilon() {return def_epsilon;}

    ///Set the default epsilon value.

    static void defaultEpsilon(Value e) {def_epsilon=e;}

    ///\name Comparisons

    ///See class Tolerance for more details.

    ///@{

    ///Returns \c true if \c a is \e surely strictly less than \c b

    bool less(Value a,Value b) const {return a+_epsilon<b;}

    ///Returns \c true if \c a is \e surely different from \c b

    bool different(Value a,Value b) const { return less(a,b)||less(b,a); }

    ///Returns \c true if \c a is \e surely positive

    bool positive(Value a) const { return _epsilon<a; }

    ///Returns \c true if \c a is \e surely negative

    bool negative(Value a) const { return -_epsilon>a; }

    ///Returns \c true if \c a is \e surely non-zero

    ///@}

    ///Returns zero

    static Value zero() {return 0;}

  };

  ///Integer specialization of \ref Tolerance.

  ///Integer specialization of \ref Tolerance.

  ///\sa Tolerance

  template<>

  class Tolerance<int>

  {

  public:

    ///\e

    typedef int Value;

    ///\name Comparisons

    ///See \ref Tolerance for more details.

    ///@{

    ///Returns \c true if \c a is \e surely strictly less than \c b

    static bool less(Value a,Value b) { return a<b;}

    ///Returns \c true if \c a is \e surely different from \c b

    static bool different(Value a,Value b) { return a!=b; }

    ///Returns \c true if \c a is \e surely positive

    static bool positive(Value a) { return 0<a; }

    ///Returns \c true if \c a is \e surely negative

    static bool negative(Value a) { return 0>a; }

    ///Returns \c true if \c a is \e surely non-zero

    ///@}

    ///Returns zero

    static Value zero() {return 0;}

  };

  ///Unsigned integer specialization of \ref Tolerance.

  ///Unsigned integer specialization of \ref Tolerance.

  ///\sa Tolerance

  template<>

  class Tolerance<unsigned int>

  {

  public:

    ///\e

    typedef unsigned int Value;

    ///\name Comparisons

    ///See \ref Tolerance for more details.

    ///@{

    ///Returns \c true if \c a is \e surely strictly less than \c b

    static bool less(Value a,Value b) { return a<b;}

    ///Returns \c true if \c a is \e surely different from \c b

    static bool different(Value a,Value b) { return a!=b; }

    ///Returns \c true if \c a is \e surely positive

    static bool positive(Value a) { return 0<a; }

    ///Returns \c true if \c a is \e surely negative

    static bool negative(Value) { return false; }

    ///Returns \c true if \c a is \e surely non-zero

    ///@}

    ///Returns zero

    static Value zero() {return 0;}

  };

  ///Long integer specialization of \ref Tolerance.

  ///Long integer specialization of \ref Tolerance.

  ///\sa Tolerance

  template<>

  class Tolerance<long int>

  {

  public:

    ///\e

    typedef long int Value;

    ///\name Comparisons

    ///See \ref Tolerance for more details.

    ///@{

    ///Returns \c true if \c a is \e surely strictly less than \c b

    static bool less(Value a,Value b) { return a<b;}

    ///Returns \c true if \c a is \e surely different from \c b

    static bool different(Value a,Value b) { return a!=b; }

    ///Returns \c true if \c a is \e surely positive

    static bool positive(Value a) { return 0<a; }

    ///Returns \c true if \c a is \e surely negative

    static bool negative(Value a) { return 0>a; }

    ///Returns \c true if \c a is \e surely non-zero

    ///@}

    ///Returns zero

    static Value zero() {return 0;}

  };

  ///Unsigned long integer specialization of \ref Tolerance.

  ///Unsigned long integer specialization of \ref Tolerance.

  ///\sa Tolerance

  template<>

  class Tolerance<unsigned long int>

  {

  public:

    ///\e

    typedef unsigned long int Value;

    ///\name Comparisons

    ///See \ref Tolerance for more details.

    ///@{

    ///Returns \c true if \c a is \e surely strictly less than \c b

    static bool less(Value a,Value b) { return a<b;}

    ///Returns \c true if \c a is \e surely different from \c b

    static bool different(Value a,Value b) { return a!=b; }

    ///Returns \c true if \c a is \e surely positive

    static bool positive(Value a) { return 0<a; }

    ///Returns \c true if \c a is \e surely negative

    static bool negative(Value) { return false; }

    ///Returns \c true if \c a is \e surely non-zero

    ///@}

    ///Returns zero

    static Value zero() {return 0;}

  };

#if defined __GNUC__ && !defined __STRICT_ANSI__

  ///Long long integer specialization of \ref Tolerance.

  ///Long long integer specialization of \ref Tolerance.

  ///\warning This class (more exactly, type <tt>long long</tt>)

  ///is not ansi compatible.

  ///\sa Tolerance

  template<>

  class Tolerance<long long int>

  {

  public:

    ///\e

    typedef long long int Value;

    ///\name Comparisons

    ///See \ref Tolerance for more details.

    ///@{

    ///Returns \c true if \c a is \e surely strictly less than \c b

    static bool less(Value a,Value b) { return a<b;}

    ///Returns \c true if \c a is \e surely different from \c b

    static bool different(Value a,Value b) { return a!=b; }

    ///Returns \c true if \c a is \e surely positive

    static bool positive(Value a) { return 0<a; }

    ///Returns \c true if \c a is \e surely negative

    static bool negative(Value a) { return 0>a; }

    ///Returns \c true if \c a is \e surely non-zero

    ///@}

    ///Returns zero

    static Value zero() {return 0;}

  };

  ///Unsigned long long integer specialization of \ref Tolerance.

  ///Unsigned long long integer specialization of \ref Tolerance.

  ///\warning This class (more exactly, type <tt>unsigned long long</tt>)

  ///is not ansi compatible.

  ///\sa Tolerance

  template<>

  class Tolerance<unsigned long long int>

  {

  public:

    ///\e

    typedef unsigned long long int Value;

    ///\name Comparisons

    ///See \ref Tolerance for more details.

    ///@{

    ///Returns \c true if \c a is \e surely strictly less than \c b

    static bool less(Value a,Value b) { return a<b;}

    ///Returns \c true if \c a is \e surely different from \c b

    static bool different(Value a,Value b) { return a!=b; }

    ///Returns \c true if \c a is \e surely positive

    static bool positive(Value a) { return 0<a; }

    ///Returns \c true if \c a is \e surely negative

    static bool negative(Value) { return false; }

    ///Returns \c true if \c a is \e surely non-zero

    ///@}

    ///Returns zero

    static Value zero() {return 0;}

  };

#endif

  /// @}

} //namespace lemon

#endif //LEMON_TOLERANCE_H