comparison sample.cpp @ 278:d9dba57becd4 sampling

Use the GNU Scientific Library implementation of MT19937 for generating
random numbers.  Answer is sufficiently similar to that from random() on
Linux that I believe I've translated things right.

 #include "audioDB.h"
 
 #include <gsl/gsl_sf.h>
+#include <gsl/gsl_rng.h>
 
 static
 double yfun(double d) {
   return gsl_sf_log(d) - gsl_sf_psi(d);
 }

   }
 
   return c;
 }
 
-unsigned audioDB::random_track(unsigned *propTable, unsigned total) {
+unsigned audioDB::random_track(unsigned *propTable, unsigned total, gsl_rng *rng) {
   /* FIXME: make this O(1) by using the alias-rejection method, or
      some other sensible method of sampling from a discrete
      distribution. */
-  /* FIXME: use a real random number generator, not random() */
-  double thing = random() / (double) RAND_MAX;
+  double thing = gsl_rng_uniform(rng);
   unsigned sofar = 0;
   for (unsigned int i = 0; i < dbH->numFiles; i++) {
     sofar += propTable[i];
     if (thing < ((double) sofar / (double) total)) {
       return i;

 }
 
 void audioDB::sample(const char *dbName) {
   initTables(dbName, 0);
 
+  gsl_rng *rng = gsl_rng_alloc(gsl_rng_mt19937);
+
   /* FIXME: in Real Life we'll want to initialize the RNG using
      /dev/random or the current time or something, like this:
 
      unsigned int seed;
      int fd = open("/dev/urandom", O_RDONLY);
      read(fd, &seed, 4);
      
-     srandom(seed);
-  */ 
+     gsl_rng_set(rng, seed);
+  */
 
   // build track offset table (FIXME: cut'n'pasted from query.cpp)
   off_t *trackOffsetTable = new off_t[dbH->numFiles];
   unsigned cumTrack=0;
   for(unsigned int k = 0; k < dbH->numFiles; k++){

 
   double sumdist = 0;
   double sumlogdist = 0;
 
   for (unsigned int i = 0; i < nsamples;) {
-    unsigned track1 = random_track(propTable, total);
-    unsigned track2 = random_track(propTable, total);
+    unsigned track1 = random_track(propTable, total, rng);
+    unsigned track2 = random_track(propTable, total, rng);
 
     if(track1 == track2)
       continue;
 
-    /* FIXME: this uses lower-order bits, which is OK on Linux but not
-       necessarily elsewhere.  Again, use a real random number
-       generator */
-    unsigned i1 = random() % propTable[track1];
-    unsigned i2 = random() % propTable[track2];
+    unsigned i1 = gsl_rng_uniform_int(rng, propTable[track1]);
+    unsigned i2 = gsl_rng_uniform_int(rng, propTable[track2]);
 
     VERB_LOG(1, "%d %d, %d %d | ", track1, i1, track2, i2);
 
     /* FIXME: this seeking, reading and distance calculation should
        share more code with the query loop */