random_util.c

关于遗传算法的一些见地。特别是关于简单遗传程序设计的实现。

  last updated:	30/04/01
 **********************************************************************/

double random_cauchy(void)
  {
  return tan(random_double_range(-PI/2,PI/2));
  }


/**********************************************************************
  random_exponential()
  synopsis:	Random number with an exponential distribution, mean
		of 1.0.
  parameters:	none
  return:	double	Random value.
  last updated:	30/04/01
 **********************************************************************/

double random_exponential(void)
  {
  return -log(random_unit_uniform());
  }


/**********************************************************************
  random_int_permutation()
  synopsis:	Randomize an array of integers.
  parameters:	int	size
		int	*iarray		Source array.
		int	*oarray		Destination array.
  return:	none
  last updated:	14 Mar 2002
 **********************************************************************/

void random_int_permutation(const int size, int *iarray, int *oarray)
  {
  int		i,j=0;		/* Loop variables over arrays. */
  int		pos;		/* Randomly selected index. */
  
  if (!iarray || !oarray) die("NULL pointer to int array passed.");

  for (i=size-1; i>0; i--)
    {
    pos = random_int(i);
    oarray[j++] = iarray[pos];
    iarray[pos] = iarray[i];
    }

  oarray[j] = iarray[0];

  return;
  }


/**********************************************************************
  random_diagnostics()
  synopsis:	Diagnostics.
  parameters:	none
  return:	none
  last updated:	13 Mar 2002
 **********************************************************************/

void random_diagnostics(void)
  {
  int	i;	/* Loop over PRNG array. */

  printf("=== PRNG routines diagnostic information =====================\n");
  printf("Version:                   %s\n", GA_VERSION_STRING);
  printf("Build date:                %s\n", GA_BUILD_DATE_STRING);
  printf("Compilation machine characteristics:\n%s\n", GA_UNAME_STRING);
  printf("--------------------------------------------------------------\n");
  printf("RANDOM_DEBUG:              %d\n", RANDOM_DEBUG);
  printf("RANDOM_RAND_MAX:           %u\n", RANDOM_RAND_MAX);
  printf("RANDOM_NUM_STATE_VALS:     %d\n", RANDOM_NUM_STATE_VALS);
#if HAVE_SLANG==1
    printf("HAVE_SLANG:                TRUE\n");
#else
    printf("HAVE_SLANG:                FALSE\n");
#endif
  printf("--------------------------------------------------------------\n");
  printf("structure                  sizeof\n");
  printf("random_state:              %lu\n", (unsigned long) sizeof(random_state));
  printf("--------------------------------------------------------------\n");

  if (is_initialised)
    {
    printf("Current state\n");
    printf("j: %d k: %d x: %d v[%d]:\n",
           current_state.j, current_state.k, current_state.x, RANDOM_NUM_STATE_VALS);
    for (i=0; i<RANDOM_NUM_STATE_VALS; i++) printf("%u ", current_state.v[i]);
    printf("\n");
    }
  else
    {
    printf("Not initialised.\n");
    }

  printf("==============================================================\n");

  return;
  }


/**********************************************************************
  random_test()
  synopsis:	Testing.
  parameters:
  return:	TRUE if all tests successful.
  last updated:	30/12/00
 **********************************************************************/

#define NUM_BINS 200
#define NUM_SAMPLES 1000000
#define NUM_CHISQ 20

boolean random_test(void)
  {
  unsigned int	i, j, k;		/* Loop variables. */
  double	r;			/* Pseudo-random number. */
  long		bins[NUM_BINS];		/* Bin. */
  double	sum=0,sumsq=0;		/* Stats. */
  int		numtrue=0, numfalse=0;	/* Count booleans. */
  unsigned int	rchi=100;		/* Number of bins in chisq test. */
  unsigned int	nchi=1000;		/* Number of samples in chisq test. */
  double	chisq;			/* Chisq error. */
  double	elimit = 2*sqrt((double)rchi);	/* Chisq error limit. */
  double	nchi_rchi = (double)nchi / (double)rchi;	/* Speed calculation. */
  FILE		*rfile=NULL;		/* Handle for file of random integers. */

  random_init();

  printf("Testing random numbers.\n");

/*
 * Uniform Distribution.
 */
  printf("Uniform distribution.  Mean should be about 0.5.\n");

  for (i=0;i<NUM_BINS;i++) bins[i] = 0;

  for (i=0;i<NUM_SAMPLES;i++)
    {
    r = random_unit_uniform();
    if (r >= 0.0 && r < 1.0)
      {
      bins[(int)(r*NUM_BINS)]++;
      sum += r;
      sumsq += SQU(r);
      }
    else
      {
      printf("Number generated out of range 0.0<=r<1.0.\n");
      }
    }
  printf("Mean = %f\n", sum / NUM_SAMPLES);
  printf("Standard deviation = %f\n", (sumsq - sum*sum/NUM_SAMPLES)/NUM_SAMPLES);

  for (i=0;i<NUM_BINS;i++)
    printf("%5.3f %ld\n", i/(double)NUM_BINS, bins[i]);

/*
 * Gaussian Distribution.
 */
  printf("Gaussian distribution.  Mean should be about 0.45.  Standard deviation should be about 0.05.\n");

  sum=0;
  sumsq=0;

  for (i=0;i<NUM_BINS;i++) bins[i] = 0;

  for (i=0;i<NUM_SAMPLES;i++)
    {
    r = random_gaussian(0.45,0.05);
    if (r >= 0.0 && r < 1.0)
      {
      bins[(int)(r*NUM_BINS)]++;
      sum += r;
      sumsq += SQU(r);
      }
    else
      {
      printf("Number generated out of range 0.0<=r<1.0.\n");
      }
    }
  printf("Mean = %f\n", sum / NUM_SAMPLES);
  printf("Standard deviation = %f\n", (sumsq - sum*sum/NUM_SAMPLES)/NUM_SAMPLES);

  for (i=0;i<NUM_BINS;i++)
    printf("%5.3f %ld\n", i/(double)NUM_BINS, bins[i]);

/*
 * Unit Gaussian Distribution.
 */
  printf("Gaussian distribution.  Mean should be about 0.0.  Standard deviation should be about 1.0.\n");

  sum=0;
  sumsq=0;

  for (i=0;i<NUM_BINS;i++) bins[i] = 0;

  for (i=0;i<NUM_SAMPLES;i++)
    {
    r = random_unit_gaussian();
    if (r >= -5.0 && r < 5.0)
      {
      bins[(int)((r+5.0)*NUM_BINS)/10]++;
      sum += r;
      sumsq += SQU(r);
      }
    else
      {
      printf("Number generated out of range -5.0<=r<5.0.\n");
      }
    }
  printf("Mean = %f\n", sum / NUM_SAMPLES);
  printf("Standard deviation = %f\n", (sumsq - sum*sum/NUM_SAMPLES)/NUM_SAMPLES);

  for (i=0;i<NUM_BINS;i++)
    printf("%5.3f %ld\n", -5.0+10*i/(double)NUM_BINS, bins[i]);

/*
 * Random Boolean.
 */
  printf("Random Booleans.  Two counts should be approximately equal.\n");

  for (i=0;i<NUM_SAMPLES;i++)
    {
    if ( random_boolean() )
      numtrue++;
    else
      numfalse++;
    }
  printf("TRUE/FALSE = %d/%d\n", numtrue, numfalse);

/*
 * Random int.
 */
  printf("Random Integers.  The distribution should be approximately uniform.\n");

  for (i=0;i<NUM_BINS;i++) bins[i] = 0;

  for (i=0;i<NUM_SAMPLES;i++)
    bins[random_int(NUM_BINS)]++;

  for (i=0;i<NUM_BINS;i++)
    printf("%u %ld\n", i, bins[i]);

/*
 * Chi squared test.  This is the standard basic test for randomness of a PRNG.
 * We would expect any moethod to fail about about one out of ten runs.
 * The error is r*t/N - N and should be within 2*sqrt(r) of r.
 */

  printf("Chi Squared Test of Random Integers.  We would expect a couple of failures.\n");

  if (rchi>NUM_BINS) die("Internal error: too few bins.");

  for (j=0;j<NUM_CHISQ;j++)
    {
    printf("Run %u. chisq should be within %f of %u.\n", j, elimit, rchi);
    for(k=0; k<10; k++)
      {
      memset(bins, 0, rchi*sizeof(long));
      chisq = 0.0;

      for(i=0; i<nchi; i++)
        bins[random_int(rchi)]++;

      for(i=0; i<rchi; i++)
        chisq += SQU((double)bins[i] - nchi_rchi);

      chisq /= nchi_rchi;

      printf("chisq = %f - %s\n", chisq, fabs(chisq - rchi)>elimit?"FAILED":"PASSED");
      }
    }

  printf("Creating file (\"randtest.dat\") of 5000 random integer numbers for external analysis.\n");

  rfile = fopen("randtest.dat", "w");

  for (i=0; i<5000; i++)
    {
    r = (double) random_rand();
    fprintf(rfile, "%f %f\n",
            /*i, r,*/
            (double)i/(double)5000, r/(double)RANDOM_RAND_MAX);
    }

  fclose(rfile);

  return TRUE;
  }


/**********************************************************************
  SLang intrinsic wrapper functions.

  FIXME: Problem with unsigned vs. signed ints.
 **********************************************************************/
#if HAVE_SLANG==1

/* These function don't need wrappers:
void random_tseed(void)
void random_init(void)
boolean random_is_init(void)
boolean random_boolean(void)
double random_unit_uniform(void)
double random_unit_gaussian(void)
void random_diagnostics(void)
boolean random_test(void)
*/

int random_rand_wrapper(void)
  {
  return (int) random_rand();
  } 


void random_seed_wrapper(int *seed)
  { 
  random_seed((unsigned int) *seed);
  return;
  } 


#if 0
THREAD_LOCK_DEFINE_STATIC(state_table_lock);
TableStruct     *state_table=NULL;        /* Table for state handles. */

int random_get_state_wrapper(void)
  {
  int	stateid;	/* State handle. */

  THREAD_LOCK(state_table_lock);
  if (!state_table) state_table=table_new();

  stateid = table_add(state_table, (vpointer) current_state);	/* FIXME: Need to copy state. */
  THREAD_UNLOCK(state_table_lock);

  return stateid;
  }


void random_set_state_wrapper(int stateid)
  {
  random_state	state;	/* State to restore. */

  THREAD_LOCK(state_table_lock);
  state = table_get_data(state_table, stateid);
  THREAD_UNLOCK(state_table_lock);

  if (!state) die("Unmatched state handle.");

  random_set_state(state);

  return;
  }
#endif


boolean random_boolean_prob_wrapper(double *prob)
  {
  return random_boolean_prob(*prob);
  }


int random_int_wrapper(int *max)
  {
  return (int) random_int((unsigned int) *max);
  }


int random_int_range_wrapper(int *min, int *max)
  {
  return random_int_range(*min, *max);
  }


double random_double_wrapper(double *max)
  {
  return random_double(*max);
  }


double random_double_range_wrapper(double *min, double *max)
  {
  return random_double_range(*min, *max);
  }


double random_gaussian_wrapper(double *mean, double *stddev)
  {
  return random_gaussian(*mean, *stddev);
  }

#endif