/*
 *  Tony Pasqualoni / Sept. 20, 2006
 *
 *  GNU Scientific Library random number generator test:
 *  Measures speed and runs Diehard tests on four GSL random number generators.
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <math.h>
#include <assert.h>
#include <time.h>

#include <gsl/gsl_rng.h>

/* 
   Compilation:
   gcc -O2 -c [filename] ; gcc -O2 [filename] -lgsl -lgslcblas -lm
 
   Using the GSL library may require the following:
   LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
   export LD_LIBRARY_PATH

   see http://www.gnu.org/software/gsl/manual/html_node/Shared-Libraries.html#Shared-Libraries
      for more information.
*/

#define MAX_STRING_SIZE 100
#define MEGA 1000000

unsigned int num_amt;

// function definitions:
void speed_test (const gsl_rng_type * gsl_type, unsigned int num_amt, unsigned int runs );
void diehard_test (const gsl_rng_type * gsl_type, unsigned int num_amt, unsigned int runs );

int main (void) {

   const gsl_rng_type **gsl_type0;
   const gsl_rng_type **gsl_type;

   // use random() for seeding RNG with random values:
   srandom(0);

   // amount of random numbers to generate for tests:
   num_amt = 10 * MEGA;

   speed_test (gsl_rng_taus, num_amt, 10);
   speed_test (gsl_rng_gfsr4, num_amt, 10);
   speed_test (gsl_rng_mt19937, num_amt, 10);
   speed_test (gsl_rng_ranlxd1, num_amt, 10);

   diehard_test (gsl_rng_taus, num_amt, 10);
   diehard_test (gsl_rng_gfsr4, num_amt, 10);
   diehard_test (gsl_rng_mt19937, num_amt, 10);
   diehard_test (gsl_rng_ranlxd1, num_amt, 10);

   exit(0);
   
   // test all generators:
   gsl_type0 = gsl_rng_types_setup();
   for (gsl_type = gsl_type0; *gsl_type != 0; gsl_type++) {
      speed_test (*gsl_type, num_amt, 1);
      diehard_test (*gsl_type, num_amt, 1);
   }

} // main

void speed_test (const gsl_rng_type * gsl_type, unsigned int num_amt, unsigned int runs ) {
// Measure speed of RNG using clock() function.

   gsl_rng * rng;
      
   int i;
   unsigned int int_32;
   clock_t start, stop;
   unsigned int seed;
   unsigned int run;
   double total_time;

   printf ("\nGenerator: %s\n\n", gsl_type->name);
   printf ("Run:   Time to generate %u million integers (seconds):\n",num_amt/MEGA);

   total_time = 0;
   for (run = 0; run < runs; run++) {
      printf("%-3d ",run);

      seed = 0;

      start = clock();
      rng = gsl_rng_alloc (gsl_type);    // initialize RNG
      gsl_rng_set (rng,seed);            // set RNG seed
      for (i = 0; i < num_amt; i++) 
         int_32 = gsl_rng_get (rng);     // generate a 32-bit unsigned integer
      stop = clock();

      total_time += ((double) (stop - start)) / CLOCKS_PER_SEC;
      printf("   %-2.6f\n",((double) (stop - start)) / CLOCKS_PER_SEC );

   } // for run

   printf("\nMean time for %s to generate %u million integers: %2.6f seconds.\n",
      gsl_type->name,num_amt/MEGA,total_time/runs);

} // void speed test

void diehard_test (const gsl_rng_type * gsl_type, unsigned int num_amt, unsigned int runs ) {
// Run Diehard tests on integers generated by RNG. 
// This function uses a Python script, 'parse-diehard.py', to parse Diehard results file.

   int i;

   unsigned int * gsl_ints;
   unsigned int run;
   
   gsl_rng * rng;

   FILE * fp;
   char fname[MAX_STRING_SIZE];
   char sys_command[MAX_STRING_SIZE];
   unsigned int seed;

   printf("\nDiehard tests for %s, %u runs:\n\n",gsl_type->name,runs);
   printf("Run:           Seed:\n");

   // back up current results file, if any:
   system ("cat 2>&1 ./result.txt >> ./result-backup.txt ; rm 2>&1 ./result.txt");

   gsl_ints = (unsigned int *) malloc (num_amt * sizeof(unsigned int));
   if (gsl_ints == NULL) {
      printf ("Malloc error.\n");
      exit(1);
   }

   // initialize RNG:
   rng = gsl_rng_alloc (gsl_type); 

   for (run = 0; run < runs; run++) {
   
      // set rng seed:
      seed = random();   
      gsl_rng_set (rng,seed);
      // generate random integers:
      for (i = 0; i < num_amt; i++) 
         gsl_ints[i] = gsl_rng_get (rng);

      // save random number array: 
      sprintf(fname,"gsl-%s-seed-%u.bit",gsl_type -> name,seed);
      printf("%3u    %12u   ",run + 1,seed);
      fp = fopen (fname,"wb");
      assert (fwrite (gsl_ints, sizeof (unsigned int), num_amt, fp));
      fclose(fp);

      // run Diehard:
      printf("running Diehard tests ... ");
      fflush(stdout);
      sprintf(sys_command,"diehard -t %s > /dev/null",fname);
      system(sys_command);

      // remove random number file:
      sprintf(sys_command,"rm %s",fname);
      system(sys_command);
      
      system("parse-diehard.py | grep -i failed");

   } 
   fflush(stdout);
   printf("\n");

   // free memory associated with rng:
   gsl_rng_free (rng);
 
   // free random number array:
   free (gsl_ints);

} // void diehard_test