source: trunk/packages/optimizer/src/pgapack/pgapack/source/random.c @ 5673

/*
 *  
 *  ********************************************************************* 
 *  (C) COPYRIGHT 1995 UNIVERSITY OF CHICAGO 
 *  *********************************************************************
 *  
 *  This software was authored by
 *  
 *  D. Levine
 *  Mathematics and Computer Science Division Argonne National Laboratory
 *  Argonne IL 60439
 *  levine@mcs.anl.gov
 *  (708) 252-6735
 *  (708) 252-5986 (FAX)
 *  
 *  with programming assistance of participants in Argonne National 
 *  Laboratory's SERS program.
 *  
 *  This program contains material protectable under copyright laws of the 
 *  United States.  Permission is hereby granted to use it, reproduce it, 
 *  to translate it into another language, and to redistribute it to 
 *  others at no charge except a fee for transferring a copy, provided 
 *  that you conspicuously and appropriately publish on each copy the 
 *  University of Chicago's copyright notice, and the disclaimer of 
 *  warranty and Government license included below.  Further, permission 
 *  is hereby granted, subject to the same provisions, to modify a copy or 
 *  copies or any portion of it, and to distribute to others at no charge 
 *  materials containing or derived from the material.
 *  
 *  The developers of the software ask that you acknowledge its use in any 
 *  document referencing work based on the  program, such as published 
 *  research.  Also, they ask that you supply to Argonne National 
 *  Laboratory a copy of any published research referencing work based on 
 *  the software.
 *  
 *  Any entity desiring permission for further use must contact:
 *  
 *  J. Gleeson
 *  Industrial Technology Development Center Argonne National Laboratory
 *  Argonne IL 60439
 *  gleesonj@smtplink.eid.anl.gov
 *  (708) 252-6055
 *  
 *  ******************************************************************** 
 *  DISCLAIMER
 *  
 *  THIS PROGRAM WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED BY AN AGENCY 
 *  OF THE UNITED STATES GOVERNMENT.  NEITHER THE UNIVERSITY OF CHICAGO, 
 *  THE UNITED STATES GOVERNMENT NOR ANY OF THEIR EMPLOYEES MAKE ANY 
 *  WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR 
 *  RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR USEFULNESS OF ANY 
 *  INFORMATION OR PROCESS DISCLOSED, OR REPRESENTS THAT ITS USE WOULD NOT 
 *  INFRINGE PRIVATELY OWNED RIGHTS.
 *  
 *  ********************************************************************** 
 *  GOVERNMENT LICENSE
 *  
 *  The Government is granted for itself and others acting on its behalf a 
 *  paid-up, non-exclusive, irrevocable worldwide license in this computer 
 *  software to reproduce, prepare derivative works, and perform publicly 
 *  and display publicly.
 */

/*****************************************************************************
*     FILE: random.c: This file contains routines to generate randomness.
*
*     Authors: David M. Levine, Philip L. Hallstrom, David M. Noelle,
*              Brian P. Walenz
*****************************************************************************/

#include "pgapack.h"

/*U****************************************************************************
   PGARandomFlip - flip a biased coin and return PGA_TRUE if the coin is
   a "winner."  Otherwise, return PGA_FALSE.

   Category: Utility

   Inputs:
      ctx - context variable
      p   - biased probability (.5 is a fair coin)

   Outputs:
      PGA_TRUE or PGA_FALSE

   Example:
      To return PGA_TRUE approximately seventy percent of the time, use

      PGAContext *ctx;
      int p;
      :
      PGARandomFlip(ctx, 0.7)

****************************************************************************U*/
int PGARandomFlip ( PGAContext *ctx, double p )
{
    PGADebugEntered("PGARandomFlip");

    PGADebugExited("PGARandomFlip");

    return( (PGARandom01(ctx, 0) < p) ? PGA_TRUE : PGA_FALSE);
}


/*U****************************************************************************
   PGARandomInterval - returns a uniform random number on the specified
   interval

   Category: Utility

   Inputs:
      ctx - context variable
      start - starting (integer) value of the interval
      end   - ending   (integer) value of the interval

   Outputs:
      A uniformly distributed random number in the interval [start, end].

   Example:
      Generate a value uniformly random from the interval [0,99]

      PGAContext *ctx;
      :
      PGARandomInterval(ctx, 0, 99);

****************************************************************************U*/
int PGARandomInterval( PGAContext *ctx, int start, int end)
{
    PGADebugEntered("PGARandomInterval");
    
    PGADebugExited("PGARandomInterval");

    return( (int)floor(PGARandom01(ctx, 0) * (double)(end-start+1) ) + start );
    
/*
   The original call...

   return(ceil( (double)(end-start+1) * PGARandom01(ctx, 0) )+start-1);
*/
}


/*****************************************************************************
*  This is a C language implementation of the universal random number        *
*  generator proposed by G. Marsaglia and A. Zaman and translated from       *
*  F. James' version.                                                        *
*                                                                            *
*  F. James                                                                  *
*  A review of pseudorandom number generators                                *
*  Computer Physics Communication                                            *
*  60 (1990) 329-344                                                         *
*                                                                            *
*  G. Marsaglia, A. Zaman, W. Tseng                                          *
*  Stat Prob. Letter                                                         *
*  9 (1990) 35.                                                              *
*                                                                            *
*  G. Marsaglia, A. Zaman                                                    *
*  FSU-SCRI-87-50                                                            *
*                                                                            *
*  This algorithm is a combination of a lagged Fibonacci and arithmetic      *
*  sequence (F. James) generator with period of 2^144.  It provides 32-bit   *
*  floating point numbers in the range from zero to one.  It is claimed to   *
*  be portable and provides bit-identical results on all machines with at    *
*  least 24-bit mantissas.                                                   *
*                                                                            *
*  PGARandom01 should be initialized with a 32-bit integer seed such that    *
*  0 <= seed <= 900,000,000.  Each of these 900,000,000 values gives rise    *
*  to an independent sequence of ~ 10^30.                                    *
*                                                                            *
*  warning on use of static storage class on thread shared memory machines   *
*****************************************************************************/
/*U****************************************************************************
   PGARandom01 - generates a uniform random number on the interval [0,1)
   If the second argument is 0 it returns the next random number in the
   sequence.  Otherwise, the second argument is used as a new seed for the
   population

   Category: Utility

   Inputs:
      ctx     - context variable
      newseed - either 0 to get the next random number, or nonzero
                to reseed
   Outputs:
      A random number on the interval [0,1)

   Example:
      To get the next random number use

      PGAContext *ctx;
      double r;
      :
      r = PGARandom01(ctx,0);

****************************************************************************U*/
double PGARandom01( PGAContext *ctx, int newseed )
{

    /* initialization variables */
    int ij, kl, i, j, k, l, m, ii, jj;
    float s, t;

    /* random number variables */
    static int seed=1;      /* default seed if none specified */
    static int i96, j96;
    static float u[97], uni, c, cd, cm;


    PGADebugEntered("PGARandom01");

    /* initialization */
/*     printf("i96 = %d\tj96 = %d\n", i96, j96); */

    if ( newseed != 0 ) {

        seed = newseed % 900000000;
        ij   = seed / 30082;
        kl   = seed - 30082 * ij;
        i    = ( (ij/177) % 177 ) + 2;
        j    = (  ij      % 177 ) + 2;
        k    = ( (kl/169) % 178 ) + 1;
        l    = (  kl      % 169 );

        for ( ii=0; ii<97; ii++ ) {

            s = 0.0;
            t = 0.5;

            for ( jj=0; jj<24; jj++ ) {

                m = ( ((i*j) % 179) * k ) % 179;
                i = j;
                j = k;
                k = m;
                l = ( (53*l) + 1 ) % 169;
                if ( ( (l*m) % 64 ) >= 32 )
                    s += t;
                t *= .5;
            }

            u[ii] = s;
        }

        c   = 362436.  /16777216.;
        cd  = 7654321. /16777216.;
        cm  = 16777213./16777216.;
        i96 = 96;
        j96 = 32;
    }

    /* random number generation */
    uni = u[i96] - u[j96];
    if ( uni < 0. ) uni += 1.0;
    u[i96] = uni;
    i96--;
    if ( i96 < 0  ) i96 = 96;
    j96--;
    if ( j96 < 0  ) j96 = 96;
    c   -= cd;
    if ( c   < 0. ) c += cm;
    uni -= c;
    if ( uni < 0. ) uni += 1.0;

    PGADebugExited("PGARandom01");
    return( (double) uni);
}

/*U****************************************************************************
   PGARandomUniform - returns a uniform random number on the interval
   [start,end]

   Category: Utility

   Inputs:
      ctx - context variable
      start - starting (double) value of the interval
      end   - ending   (double) value of the interval

   Outputs:
      A random number on the interval [start,end]

   Example:
      Generate a uniform random number on the interval [-0.5, 1.5]

      PGAContext *ctx;
      double r;
      :
      r = PGARandomUniform(ctx, -0.5, 1.5);

****************************************************************************U*/
double PGARandomUniform( PGAContext *ctx, double start, double end)
{
    double val, r;

    PGADebugEntered("PGARandomUniform");

    r = PGARandom01(ctx, 0);
    val = (end-start) * r + start;

    PGADebugExited("PGARandomUniform");

    return(val);
}


/*U****************************************************************************
   PGARandomGaussian - returns an approximation to a Gaussian random number

   Category: Utility

   Inputs:
       mean  - the mean of the Gaussian distribution
       sigma - the standard deviation of the Gaussian distribution

   Outputs:
      A random number selected from a Gaussian distribution with given
      mean and standard deviation

   Example:
      To generate a Gaussian random number with mean 0.0 and standard
      deviation 1.0 use

      PGAContext *ctx;
      :
      r = PGARandomGaussian(ctx, 0.0, 1.0);

****************************************************************************U*/
double PGARandomGaussian( PGAContext *ctx, double mean, double sigma)
{
    int i;
    double sum = 0.;

    PGADebugEntered("PGARandomGaussian");

    for (i=11;i>=0; i--)
        sum += PGARandom01(ctx, 0);

    PGADebugExited("PGARandomGaussian");

    return ( (sum-6.0) * sigma + mean );
}

/*U***************************************************************************
   PGAGetRandomSeed - returns the integer to seed random numbers with

   Category: Utility

   Inputs:
      ctx - context variable

   Outputs:
      The seed for the random number generator

   Example:
      PGAContext *ctx;
      int seed;
      :
      seed = PGAGetRandomSeed(ctx);

***************************************************************************U*/
int PGAGetRandomSeed(PGAContext *ctx)
{
    PGADebugEntered("PGAGetRandomSeed");

    PGADebugExited("PGAGetRandomSeed");

    return(ctx->init.RandomSeed);
}

/*U****************************************************************************
   PGASetRandomSeed - set a seed for the random number generator.  The
   default is to use a random seed.  Specifying a seed exlicitly allows
   for reproducibility of runs.

   Category: Utility

   Inputs:
      ctx  - context variable
      seed - seed  for the random number generator

   Outputs:
      None

   Example:
      PGAContext *ctx;
      :
      PGASetRandomSeed(ctx,1);

****************************************************************************U*/
void PGASetRandomSeed(PGAContext *ctx, int seed)
{
#define MAX_PROCESSORS 2048

    PGADebugEntered("PGASetRandomSeed");
    PGAFailIfSetUp("PGASetRandomSeed");

    if ((seed < 1) || (seed + MAX_PROCESSORS > 900000000))
        PGAError ( ctx, "PGASetRandomSeed: Invalid value of seed:",
                  PGA_FATAL, PGA_INT, (void *) &seed);
    else
        ctx->init.RandomSeed = seed;
    
    PGADebugExited("PGASetRandomSeed");
}

/*****************************************************************************
 Category: Generation/Utility
 
 Inputs:
                ctx - context variable
                randreplaceprop - proportion of individuals to be replaced per generation randomly, post initialization
                
 Outputs:
                None
                
 Example:
        PGAContext *ctx;
        .
        .
        PGASetRandReplProp(ctx, 0.05);          
 ***************************************************************************/
void PGASetRandReplProp(PGAContext *ctx, double randreplprop){
        PGADebugEntered("PGASetRandReplProp");
        PGAFailIfSetUp("PGASetRandReplProp");
        if(randreplprop < 0 || randreplprop > 1){
                PGAError ( ctx, "PGASetRandReplProp: Invalid value set:",
                  PGA_FATAL, PGA_DOUBLE, (void *) &randreplprop);
        }else{
                ctx->ga.randreplprop = randreplprop;
        }
        PGADebugExited("PGASetRandReplProp");
}

/*****************************************************************************
 Category: Generation/Utility
 
 Inputs:
                ctx - context variable
                
 Outputs:
                randreplprop - proportion of randomly generated individuals per population, post initialization
                
 Example:
        PGAContext *ctx;
        double randreplprop;
        .
        .
        randreplprop = PGAGetRandReplProp(ctx);         
 ***************************************************************************/
double PGAGetRandReplProp(PGAContext *ctx){
        PGADebugEntered("PGAGetRandReplProp");
        PGADebugExited("PGAGetRandReplProp");
        return ctx->ga.randreplprop;
}