[Commit-gnuradio] [gnuradio] 03/13: add boost.random as random number generator

[git]

This is an automated email from the git hooks/post-receive script.

jcorgan pushed a commit to branch master
in repository gnuradio.

commit e172d55e2bfe63e5b99e7d432f9386c37bff8c84
Author: Stefan <address@hidden>
Date:   Tue Sep 1 15:09:17 2015 +0200

    add boost.random as random number generator
---
 gnuradio-runtime/include/gnuradio/random.h       |  21 +++--
 gnuradio-runtime/lib/math/random.cc              | 106 +++++++++--------------
 gnuradio-runtime/python/gnuradio/gr/qa_random.py |  19 ++++
 3 files changed, 73 insertions(+), 73 deletions(-)

diff --git a/gnuradio-runtime/include/gnuradio/random.h 
b/gnuradio-runtime/include/gnuradio/random.h
index 454c88e..f15b12a 100644
--- a/gnuradio-runtime/include/gnuradio/random.h
+++ b/gnuradio-runtime/include/gnuradio/random.h
@@ -37,6 +37,8 @@ static const int RANDOM_MAX = 2147483647; // 2^31-1
 #endif /* RANDOM_MAX */
 
 #include <stdlib.h>
+#include <boost/random.hpp>
+#include <ctime>
 
 namespace gr {
 
@@ -47,17 +49,22 @@ namespace gr {
   class GR_RUNTIME_API random
   {
   protected:
-    static const int NTAB = 32;
     long d_seed;
-    long d_iy;
-    long d_iv[NTAB];
-    int        d_iset;
-    float d_gset;
+    bool d_gauss_stored;
+    float d_gauss_value;
+
+    boost::mt19937 *d_rng; // mersenne twister as random number generator
+    boost::uniform_real<float> *d_uniform; // choose uniform distribution, 
default is [0,1)
+    boost::variate_generator<boost::mt19937&, boost::uniform_real<float> > 
*d_generator;
 
   public:
-    random(long seed=3021);
+    random(unsigned int seed=0);
+    ~random();
 
-    void reseed(long seed);
+    /*!
+     * \brief Change the seed for the initialized number generator
+     */
+    void reseed(unsigned int seed);
 
     /*!
      * \brief Uniform random numbers in the range [0.0, 1.0)
diff --git a/gnuradio-runtime/lib/math/random.cc 
b/gnuradio-runtime/lib/math/random.cc
index ecf70e0..76de3bd 100644
--- a/gnuradio-runtime/lib/math/random.cc
+++ b/gnuradio-runtime/lib/math/random.cc
@@ -44,90 +44,67 @@
 
 namespace gr {
 
-#define IA 16807
-#define IM 2147483647
-#define AM (1.0/IM)
-#define IQ 127773
-#define IR 2836
-#define NDIV (1+(IM-1)/NTAB)
-#define EPS 1.2e-7
-#define RNMX (1.0-EPS)
+  random::random(unsigned int seed)
+  {
+    d_gauss_stored = false; // set gasdev (gauss distributed numbers) on 
calculation state
+
+    // Setup random number generator
+    d_rng = new boost::mt19937;
+    reseed(seed); // set seed for random number generator
+    d_uniform = new boost::uniform_real<float>;
+    d_generator = new boost::variate_generator<boost::mt19937&, 
boost::uniform_real<float> > (*d_rng,*d_uniform); // create number generator in 
[0,1) from boost.random
+  }
 
-  random::random(long seed)
+  random::~random()
   {
-    reseed(seed);
+      delete d_rng;
+      delete d_uniform;
+      delete d_generator;
   }
 
+  /*
+   * Seed is initialized with time if the given seed is 0. Otherwise the seed 
is taken directly. Sets the seed for the random number generator.
+   */
   void
-  random::reseed(long seed)
+  random::reseed(unsigned int seed)
   {
-    d_seed = seed;
-    d_iy = 0;
-    for(int i = 0; i < NTAB; i++)
-      d_iv[i] = 0;
-    d_iset = 0;
-    d_gset = 0;
+    if(seed==0) d_seed = static_cast<unsigned int>(std::time(0)); // FIXME: 
add seed method correctly
+    else d_seed = seed;
+    d_rng->seed(d_seed);
   }
 
   /*
-   * This looks like it returns a uniform random deviate between 0.0 and 1.0
-   * It looks similar to code from "Numerical Recipes in C".
+   * Returns uniformly distributed numbers in [0,1) taken from boost.random 
using a Mersenne twister
    */
   float
   random::ran1()
   {
-    int j;
-    long k;
-    float temp;
-
-    if(d_seed <= 0 || !d_iy)  {
-      if(-d_seed < 1)
-        d_seed=1;
-      else
-        d_seed = -d_seed;
-      for(j=NTAB+7;j>=0;j--) {
-        k=d_seed/IQ;
-        d_seed=IA*(d_seed-k*IQ)-IR*k;
-        if(d_seed < 0)
-          d_seed += IM;
-        if(j < NTAB)
-          d_iv[j] = d_seed;
-      }
-      d_iy=d_iv[0];
-    }
-    k=(d_seed)/IQ;
-    d_seed=IA*(d_seed-k*IQ)-IR*k;
-    if(d_seed < 0)
-      d_seed += IM;
-    j=d_iy/NDIV;
-    d_iy=d_iv[j];
-    d_iv[j] = d_seed;
-    temp=AM * d_iy;
-    if(temp > RNMX)
-      temp = RNMX;
-    return temp;
+    return (*d_generator)();
   }
 
   /*
    * Returns a normally distributed deviate with zero mean and variance 1.
-   * Also looks like it's from "Numerical Recipes in C".
+   * Used is the Marsaglia polar method.
+   * Every second call a number is stored because the transformation works 
only in pairs. Otherwise half calculation is thrown away.
    */
   float
   random::gasdev()
   {
-    float fac,rsq,v1,v2;
-    d_iset = 1 - d_iset;
-    if(d_iset) {
-      do {
-        v1=2.0*ran1()-1.0;
-        v2=2.0*ran1()-1.0;
-        rsq=v1*v1+v2*v2;
-      } while(rsq >= 1.0 || rsq == 0.0);
-      fac= sqrt(-2.0*log(rsq)/rsq);
-      d_gset=v1*fac;
-      return v2*fac;
-    }
-    return d_gset;
+      if(d_gauss_stored){ // just return the stored value if available
+          d_gauss_stored = false;
+          return d_gauss_value;
+      }
+      else{ // generate a pair of gaussian distributed numbers
+          float x,y,s;
+          do{
+              x = 2.0*ran1()-1.0;
+              y = 2.0*ran1()-1.0;
+              s = x*x+y*y;
+          }while(not(s<1.0));
+          d_gauss_stored = true;
+          d_gauss_value = x*sqrt(-2.0*log(s)/s);
+          return y*sqrt(-2.0*log(s)/s);
+      }
   }
 
   float
@@ -153,9 +130,6 @@ namespace gr {
       return z;
   }
 
-  /*
-   * Complex rayleigh is really gaussian I and gaussian Q
-   */
   gr_complex
   random::rayleigh_complex()
   {
diff --git a/gnuradio-runtime/python/gnuradio/gr/qa_random.py 
b/gnuradio-runtime/python/gnuradio/gr/qa_random.py
index ee40183..ef85337 100644
--- a/gnuradio-runtime/python/gnuradio/gr/qa_random.py
+++ b/gnuradio-runtime/python/gnuradio/gr/qa_random.py
@@ -110,5 +110,24 @@ class test_random(gr_unittest.TestCase):
         for k in range(len(hist[0])):
             print hist[1][k], hist[1][k+1], hist[0][k], 
float(rayleigh.cdf(hist[1][k+1])-rayleigh.cdf(hist[1][k]))*self.num_tests
 
+    # Check seeds (init with time and seed as fix number, no assert)
+    def test_6(self):
+        print '# TEST 6'
+        rndm0 = gr.random(0); # init with time
+        rndm1 = gr.random(42); # init with fix seed
+        num = 5
+
+        print 'Some random numbers in [0,1), should change every run:'
+        for k in range(num):
+            print rndm0.ran1(),
+        print ' '
+
+        print 'Some random numbers in [0,1), should be the same every run:'
+        for k in range(num):
+            print rndm1.ran1(),
+        print '== '
+        print '0.374540120363 0.796543002129 0.950714290142 0.183434784412 
0.731993913651'
+
+
 if __name__ == '__main__':
     gr_unittest.run(test_random, "test_random.xml")