COIN-OR::LEMON - Graph Library

Changeset 209:765619b7cbb2 in lemon-1.0 for lemon/random.h


Ignore:
Timestamp:
07/13/08 20:51:02 (11 years ago)
Author:
Alpar Juttner <alpar@…>
Branch:
default
Phase:
public
Message:

Apply unify-sources.sh to the source tree

File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/random.h

    r178 r209  
    1 /* -*- C++ -*-
    2  *
    3  * This file is a part of LEMON, a generic C++ optimization library
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
    44 *
    55 * Copyright (C) 2003-2008
     
    2222 *
    2323 * See the appropriate copyright notice below.
    24  * 
     24 *
    2525 * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
    26  * All rights reserved.                         
     26 * All rights reserved.
    2727 *
    2828 * Redistribution and use in source and binary forms, with or without
     
    3737 *    documentation and/or other materials provided with the distribution.
    3838 *
    39  * 3. The names of its contributors may not be used to endorse or promote 
    40  *    products derived from this software without specific prior written 
     39 * 3. The names of its contributors may not be used to endorse or promote
     40 *    products derived from this software without specific prior written
    4141 *    permission.
    4242 *
     
    8888
    8989  namespace _random_bits {
    90    
     90
    9191    template <typename _Word, int _bits = std::numeric_limits<_Word>::digits>
    9292    struct RandomTraits {};
     
    100100      static const int length = 624;
    101101      static const int shift = 397;
    102      
     102
    103103      static const Word mul = 0x6c078965u;
    104104      static const Word arrayInit = 0x012BD6AAu;
     
    168168          0x12345u, 0x23456u, 0x34567u, 0x45678u
    169169        };
    170    
     170
    171171        initState(seedArray, seedArray + 4);
    172172      }
     
    176176        static const Word mul = RandomTraits<Word>::mul;
    177177
    178         current = state; 
     178        current = state;
    179179
    180180        Word *curr = state + length - 1;
     
    202202        num = length > end - begin ? length : end - begin;
    203203        while (num--) {
    204           curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1)) 
     204          curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1))
    205205            + *it + cnt;
    206206          ++it; ++cnt;
     
    224224          }
    225225        }
    226        
     226
    227227        state[length - 1] = Word(1) << (bits - 1);
    228228      }
    229      
     229
    230230      void copyState(const RandomCore& other) {
    231231        std::copy(other.state, other.state + length, state);
     
    242242    private:
    243243
    244  
     244
    245245      void fillState() {
    246246        static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask };
     
    248248        static const Word hiMask = RandomTraits<Word>::hiMask;
    249249
    250         current = state + length; 
     250        current = state + length;
    251251
    252252        register Word *curr = state + length - 1;
    253253        register long num;
    254      
     254
    255255        num = length - shift;
    256256        while (num--) {
     
    270270      }
    271271
    272  
     272
    273273      Word *current;
    274274      Word state[length];
    275      
    276     };
    277 
    278 
    279     template <typename Result, 
     275
     276    };
     277
     278
     279    template <typename Result,
    280280              int shift = (std::numeric_limits<Result>::digits + 1) / 2>
    281281    struct Masker {
     
    285285      }
    286286    };
    287    
     287
    288288    template <typename Result>
    289289    struct Masker<Result, 1> {
     
    293293    };
    294294
    295     template <typename Result, typename Word, 
    296               int rest = std::numeric_limits<Result>::digits, int shift = 0, 
     295    template <typename Result, typename Word,
     296              int rest = std::numeric_limits<Result>::digits, int shift = 0,
    297297              bool last = rest <= std::numeric_limits<Word>::digits>
    298298    struct IntConversion {
    299299      static const int bits = std::numeric_limits<Word>::digits;
    300    
     300
    301301      static Result convert(RandomCore<Word>& rnd) {
    302302        return static_cast<Result>(rnd() >> (bits - rest)) << shift;
    303303      }
    304      
    305     }; 
    306 
    307     template <typename Result, typename Word, int rest, int shift> 
     304
     305    };
     306
     307    template <typename Result, typename Word, int rest, int shift>
    308308    struct IntConversion<Result, Word, rest, shift, false> {
    309309      static const int bits = std::numeric_limits<Word>::digits;
    310310
    311311      static Result convert(RandomCore<Word>& rnd) {
    312         return (static_cast<Result>(rnd()) << shift) | 
     312        return (static_cast<Result>(rnd()) << shift) |
    313313          IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd);
    314314      }
     
    317317
    318318    template <typename Result, typename Word,
    319               bool one_word = (std::numeric_limits<Word>::digits < 
    320                                std::numeric_limits<Result>::digits) >
     319              bool one_word = (std::numeric_limits<Word>::digits <
     320                               std::numeric_limits<Result>::digits) >
    321321    struct Mapping {
    322322      static Result map(RandomCore<Word>& rnd, const Result& bound) {
     
    325325        Result num;
    326326        do {
    327           num = IntConversion<Result, Word>::convert(rnd) & mask; 
     327          num = IntConversion<Result, Word>::convert(rnd) & mask;
    328328        } while (num > max);
    329329        return num;
     
    351351        res *= res;
    352352        if ((exp & 1) == 1) res *= static_cast<Result>(2.0);
    353         return res; 
     353        return res;
    354354      }
    355355    };
     
    361361        res *= res;
    362362        if ((exp & 1) == 1) res *= static_cast<Result>(0.5);
    363         return res; 
     363        return res;
    364364      }
    365365    };
     
    368368    struct ShiftMultiplier<Result, 0, true> {
    369369      static const Result multiplier() {
    370         return static_cast<Result>(1.0); 
     370        return static_cast<Result>(1.0);
    371371      }
    372372    };
     
    375375    struct ShiftMultiplier<Result, -20, true> {
    376376      static const Result multiplier() {
    377         return static_cast<Result>(1.0/1048576.0); 
    378       }
    379     };
    380    
     377        return static_cast<Result>(1.0/1048576.0);
     378      }
     379    };
     380
    381381    template <typename Result>
    382382    struct ShiftMultiplier<Result, -32, true> {
    383383      static const Result multiplier() {
    384         return static_cast<Result>(1.0/424967296.0); 
     384        return static_cast<Result>(1.0/424967296.0);
    385385      }
    386386    };
     
    389389    struct ShiftMultiplier<Result, -53, true> {
    390390      static const Result multiplier() {
    391         return static_cast<Result>(1.0/9007199254740992.0); 
     391        return static_cast<Result>(1.0/9007199254740992.0);
    392392      }
    393393    };
     
    396396    struct ShiftMultiplier<Result, -64, true> {
    397397      static const Result multiplier() {
    398         return static_cast<Result>(1.0/18446744073709551616.0); 
     398        return static_cast<Result>(1.0/18446744073709551616.0);
    399399      }
    400400    };
     
    408408
    409409    template <typename Result, typename Word,
    410               int rest = std::numeric_limits<Result>::digits, int shift = 0, 
     410              int rest = std::numeric_limits<Result>::digits, int shift = 0,
    411411              bool last = rest <= std::numeric_limits<Word>::digits>
    412     struct RealConversion{ 
     412    struct RealConversion{
    413413      static const int bits = std::numeric_limits<Word>::digits;
    414414
     
    420420
    421421    template <typename Result, typename Word, int rest, int shift>
    422     struct RealConversion<Result, Word, rest, shift, false> { 
     422    struct RealConversion<Result, Word, rest, shift, false> {
    423423      static const int bits = std::numeric_limits<Word>::digits;
    424424
     
    459459      Word buffer;
    460460      int num;
    461      
     461
    462462      BoolProducer() : num(0) {}
    463463
     
    530530    // Architecture word
    531531    typedef unsigned long Word;
    532    
     532
    533533    _random_bits::RandomCore<Word> core;
    534534    _random_bits::BoolProducer<Word> bool_producer;
    535    
     535
    536536
    537537  public:
     
    555555    /// to the architecture word type.
    556556    template <typename Number>
    557     Random(Number seed) { 
     557    Random(Number seed) {
    558558      _random_bits::Initializer<Number, Word>::init(core, seed);
    559559    }
     
    565565    /// architecture word type.
    566566    template <typename Iterator>
    567     Random(Iterator begin, Iterator end) { 
     567    Random(Iterator begin, Iterator end) {
    568568      typedef typename std::iterator_traits<Iterator>::value_type Number;
    569569      _random_bits::Initializer<Number, Word>::init(core, begin, end);
     
    598598    /// converted to the architecture word type.
    599599    template <typename Number>
    600     void seed(Number seed) { 
     600    void seed(Number seed) {
    601601      _random_bits::Initializer<Number, Word>::init(core, seed);
    602602    }
     
    608608    /// architecture word type.
    609609    template <typename Iterator>
    610     void seed(Iterator begin, Iterator end) { 
     610    void seed(Iterator begin, Iterator end) {
    611611      typedef typename std::iterator_traits<Iterator>::value_type Number;
    612612      _random_bits::Initializer<Number, Word>::init(core, begin, end);
     
    626626      return false;
    627627    }
    628    
     628
    629629    /// \brief Seeding from file
    630630    ///
     
    641641    /// \return True when the seeding successes.
    642642#ifndef WIN32
    643     bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0) 
     643    bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)
    644644#else
    645     bool seedFromFile(const std::string& file = "", int offset = 0) 
     645    bool seedFromFile(const std::string& file = "", int offset = 0)
    646646#endif
    647647    {
     
    661661    /// random sequence.
    662662    /// \return Currently always true.
    663     bool seedFromTime() {       
     663    bool seedFromTime() {
    664664#ifndef WIN32
    665665      timeval tv;
     
    697697    /// It returns a random real number from the range [0, b).
    698698    template <typename Number>
    699     Number real(Number b) { 
    700       return real<Number>() * b; 
     699    Number real(Number b) {
     700      return real<Number>() * b;
    701701    }
    702702
     
    705705    /// It returns a random real number from the range [a, b).
    706706    template <typename Number>
    707     Number real(Number a, Number b) { 
    708       return real<Number>() * (b - a) + a; 
     707    Number real(Number a, Number b) {
     708      return real<Number>() * (b - a) + a;
    709709    }
    710710
     
    726726    /// It returns a random real number from the range [0, b).
    727727    template <typename Number>
    728     Number operator()(Number b) { 
    729       return real<Number>() * b; 
     728    Number operator()(Number b) {
     729      return real<Number>() * b;
    730730    }
    731731
     
    734734    /// It returns a random real number from the range [a, b).
    735735    template <typename Number>
    736     Number operator()(Number a, Number b) { 
    737       return real<Number>() * (b - a) + a; 
     736    Number operator()(Number a, Number b) {
     737      return real<Number>() * (b - a) + a;
    738738    }
    739739
     
    785785    template <typename Number>
    786786    Number integer() {
    787       static const int nb = std::numeric_limits<Number>::digits + 
     787      static const int nb = std::numeric_limits<Number>::digits +
    788788        (std::numeric_limits<Number>::is_signed ? 1 : 0);
    789789      return _random_bits::IntConversion<Number, Word, nb>::convert(core);
     
    793793      return integer<int>();
    794794    }
    795    
     795
    796796    /// \brief Returns a random bool
    797797    ///
     
    807807    ///\name Non-uniform distributions
    808808    ///
    809    
     809
    810810    ///@{
    811    
     811
    812812    /// \brief Returns a random bool
    813813    ///
     
    823823    /// transformation is used to generate a random normal distribution.
    824824    /// \todo Consider using the "ziggurat" method instead.
    825     double gauss() 
     825    double gauss()
    826826    {
    827827      double V1,V2,S;
    828828      do {
    829         V1=2*real<double>()-1;
    830         V2=2*real<double>()-1;
    831         S=V1*V1+V2*V2;
     829        V1=2*real<double>()-1;
     830        V2=2*real<double>()-1;
     831        S=V1*V1+V2*V2;
    832832      } while(S>=1);
    833833      return std::sqrt(-2*std::log(S)/S)*V1;
     
    855855
    856856    /// This function generates a gamma distribution random number.
    857     /// 
     857    ///
    858858    ///\param k shape parameter (<tt>k>0</tt> integer)
    859     double gamma(int k) 
     859    double gamma(int k)
    860860    {
    861861      double s = 0;
     
    863863      return s;
    864864    }
    865    
     865
    866866    /// Gamma distribution with given shape and scale parameter
    867867
    868868    /// This function generates a gamma distribution random number.
    869     /// 
     869    ///
    870870    ///\param k shape parameter (<tt>k>0</tt>)
    871871    ///\param theta scale parameter
     
    877877      const double v0=E/(E-delta);
    878878      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           }
     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          }
    892892      } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
    893893      return theta*(xi+gamma(int(std::floor(k))));
    894894    }
    895    
     895
    896896    /// Weibull distribution
    897897
    898898    /// This function generates a Weibull distribution random number.
    899     /// 
     899    ///
    900900    ///\param k shape parameter (<tt>k>0</tt>)
    901901    ///\param lambda scale parameter (<tt>lambda>0</tt>)
     
    904904    {
    905905      return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
    906     } 
    907      
     906    }
     907
    908908    /// Pareto distribution
    909909
    910910    /// This function generates a Pareto distribution random number.
    911     /// 
     911    ///
    912912    ///\param k shape parameter (<tt>k>0</tt>)
    913913    ///\param x_min location parameter (<tt>x_min>0</tt>)
     
    916916    {
    917917      return exponential(gamma(k,1.0/x_min))+x_min;
    918     } 
    919      
     918    }
     919
    920920    /// Poisson distribution
    921921
    922922    /// This function generates a Poisson distribution random number with
    923923    /// parameter \c lambda.
    924     /// 
     924    ///
    925925    /// The probability mass function of this distribusion is
    926926    /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]
     
    928928    /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the
    929929    /// return value.
    930    
     930
    931931    int poisson(double lambda)
    932932    {
     
    935935      double p = 1.0;
    936936      do {
    937         k++;
    938         p*=real<double>();
     937        k++;
     938        p*=real<double>();
    939939      } while (p>=l);
    940940      return k-1;
    941     } 
    942      
     941    }
     942
    943943    ///@}
    944    
     944
    945945    ///\name Two dimensional distributions
    946946    ///
    947947
    948948    ///@{
    949    
     949
    950950    /// Uniform distribution on the full unit circle
    951951
    952952    /// Uniform distribution on the full unit circle.
    953953    ///
    954     dim2::Point<double> disc() 
     954    dim2::Point<double> disc()
    955955    {
    956956      double V1,V2;
    957957      do {
    958         V1=2*real<double>()-1;
    959         V2=2*real<double>()-1;
    960        
     958        V1=2*real<double>()-1;
     959        V2=2*real<double>()-1;
     960
    961961      } while(V1*V1+V2*V2>=1);
    962962      return dim2::Point<double>(V1,V2);
     
    974974      double V1,V2,S;
    975975      do {
    976         V1=2*real<double>()-1;
    977         V2=2*real<double>()-1;
    978         S=V1*V1+V2*V2;
     976        V1=2*real<double>()-1;
     977        V2=2*real<double>()-1;
     978        S=V1*V1+V2*V2;
    979979      } while(S>=1);
    980980      double W=std::sqrt(-2*std::log(S)/S);
     
    985985    /// This function provides a turning symmetric two-dimensional distribution.
    986986    /// The x-coordinate is of conditionally exponential distribution
    987     /// with the condition that x is positive and y=0. If x is negative and 
     987    /// with the condition that x is positive and y=0. If x is negative and
    988988    /// y=0 then, -x is of exponential distribution. The same is true for the
    989989    /// y-coordinate.
    990     dim2::Point<double> exponential2() 
     990    dim2::Point<double> exponential2()
    991991    {
    992992      double V1,V2,S;
    993993      do {
    994         V1=2*real<double>()-1;
    995         V2=2*real<double>()-1;
    996         S=V1*V1+V2*V2;
     994        V1=2*real<double>()-1;
     995        V2=2*real<double>()-1;
     996        S=V1*V1+V2*V2;
    997997      } while(S>=1);
    998998      double W=-std::log(S)/S;
     
    10001000    }
    10011001
    1002     ///@}   
     1002    ///@}
    10031003  };
    10041004
Note: See TracChangeset for help on using the changeset viewer.