gen_init.c

统计模式识别算法包

      if (fit < minfitness)
         minfitness = fit ;
      }

   avgfitness /= (double) popsize ;

/*
   The second step is to apply a linear transform to these fitnesses to prevent
   extraordinarily fit individuals from dominating early on, and at the same
   time still favor the most fit later in the run when a large number of
   individuals are very fit.
   
   This transform is:  f' = tmult * f + tconst.
   
   The coefficients are chosen so that the transformed maximum fitness is
   favor_best times the transformed average, while the average after transform
   is equal to that before.  A typical value for favor_best is 2-3.
   One problem is that late in the run, when the average is close to the max,
   very small fitnesses may map negative.  In this case, map the smallest
   to zero and do the best we can for the max.

   Note that a common alternative is to use the mapping just described, and
   truncate transformed fitnesses at zero.  However, the method shown here
   is usually superior, as it preserves genetic diversity.
*/

   ftemp = maxfitness - avgfitness ;
   if (ftemp > 1.e-20) {  // Insurance: average may equal max!
      tmult = (favor_best - 1.0) * avgfitness / ftemp ;
      tconst = avgfitness * (maxfitness - favor_best * avgfitness) / ftemp ;
      }
   else {
      tmult = 1.0 ;
      tconst = 0.0 ;
      }

/*
   The 'ideal' scaling factor was just computed.  Use it to map the minimum
   fitness.  If it comes out negative, compute an alternative scaling factor
   which will map the minimum to zero and keep the average unchanged.
*/

   if (tmult * minfitness + tconst < 0.0) { // Do not allow negative fitness
      ftemp = avgfitness - minfitness ;
      if (ftemp > 1.e-20) {
         tmult = avgfitness / ftemp ;
         tconst = -minfitness * avgfitness / ftemp ;
         }
      else {
         tmult = 1.0 ;
         tconst = 0.0 ;
         }
      }

/*
   The scaling factors have been computed.  Do the scaling now.
   The truncation at zero is theoretically unnecessary, as we avoided
   negatives when we computed the scaling factor above.  However, floating
   point problems can sometimes still produce a 'negative zero'.  In deference
   to possible user modifications which DEMAND nonnegative fitnesses, it is
   good programming practice to enforce this.
*/

   avgfitness = 0.0 ;
   for (individual=0 ; individual<popsize ; individual++) {
      fit = tmult * fitness[individual] + tconst ;
      if (fit < 0.0)
         fit = 0.0 ;
      fitness[individual] = fit ;
      avgfitness += fit ;
      }

   avgfitness /= (double) popsize ;

/*
   The final step is to normalize the fitnesses by dividing each by the
   average fitness.  The effect is that then each fitness can be interpreted
   as the expected number of times it would be chosen from the population
   pool if its probability of selection were proportional to its fitness.
*/

   for (individual=0 ; individual<popsize ; individual++)
      fitness[individual] /= avgfitness ;
}

/*
--------------------------------------------------------------------------------

   fitness_to_choices - Convert the array of fitnesses (which contain
               expected frequency of selection) into the array of parent
               choices.  This will allow random selection of parents
               without replacement later, while still insuring that
               we select (to within one) the expected number of each.


--------------------------------------------------------------------------------
*/

static void fitness_to_choices (
   int popsize ,      // Length of fitness, choices vectors
   double *fitness ,  // Input array of expected selection frequencies
   int *choices       // Output array of parents
   )
{
   int individual, expected, k ;
   double rn ;

/*
   We build the choices array in two steps.  This, the first step, assigns
   parents according to the integer part of their expected frequencies.
*/

   k = 0 ;  // Will index choices array
   for (individual=0 ; individual<popsize ; individual++) {
      expected = (int) fitness[individual] ; // Assign this many now
      fitness[individual] -= expected ;      // Save fractional remainder
      while (expected--)                     // Forcibly use the int expected
         choices[k++] = individual ;         // quantity of this individual
      }

/*
   The second step is to take care of the remaining fractional expected
   frequencies.  Pass through the population, randomly selecting members
   with probability equal to their remaining fractional expectation.
   It is tempting to think that the algorithm below could loop excessively
   due to a very small fitness.  But recall that the sum of the fitnesses will
   be AT LEAST as large as the number remaining to be selected, and generally
   much more.  Thus, the ones with very small fitness (likely to cause trouble)
   will never become the only remaining possibilities.
*/

   while (k < popsize) {  // Select until choices is full
      individual = unifrand() * popsize ;// Randomly select individual
      if (fitness[individual] > 0.0) {   // Try members still having expectation
         if (fitness[individual] >= unifrand()) { // Selects with this prob
            choices[k++] = individual ;   // Bingo!  Select this individual
            fitness[individual] -= 1.0 ;  // and make it ineligable for future
            }
         }
      }
}

/*
--------------------------------------------------------------------------------

   pick_parents - Randomly select two parents from 'choices' candidate array

--------------------------------------------------------------------------------
*/

static void pick_parents (
   int *nchoices ,  // Number of choices (returned decremented by two)
   int *choices ,   // Array (nchoices long) of candidates for parent
   int *parent1 ,   // One parent returned here
   int *parent2     // and the other here
   )
{
   int k ;

   k = unifrand() * *nchoices ;        // Select position in choices array
   *parent1 = choices[k] ;             // Then return that parent
   choices[k] = choices[--*nchoices] ; // without replacement

   k = unifrand() * *nchoices ;
   *parent2 = choices[k] ;
   choices[k] = choices[--*nchoices] ;
}

/*
--------------------------------------------------------------------------------

   reproduce - Create a child from half of each of two parents.
      If first_child is true, randomly generate the crossover point.
      Else use the supplied crossover point and take other halves.
      If the chromosome size is at least 16, use two point crossover
      half of the time.  Signal this by returning a negative crosspt.
      We implement splitting the chromosome at a bit level (i.e. within
      a byte) with the 'split' parameter.  This is in the range 0-7,
      and is the bit number where the split takes place.  Like crosspt
      it is output for the first child and input for the second.

--------------------------------------------------------------------------------
*/

static void reproduce (
   char *p1 ,        // Pointer to one parent
   char *p2 ,        // and the other
   int first_child , // Is this the first of their 2 children?
   int chromsize ,   // Number of genes in chromosome
   char *child ,     // Output of a child
   int *crosspt ,    // If first_child, output of xover pt, else input it.
   int *split        // In/out of within byte splitting point
   )

{
   int i, n1, n2, n3, n4 ;
   char left, right, *pa, *pb ;

   if (first_child) {

      *split = longrand() % 8 ; // We will split boundary bytes here
      *crosspt = 1 + unifrand() * chromsize ;  // Randomly select cross pt

      if ((chromsize >= 16)  &&  (unifrand() < 0.33333)) // Two point?
         *crosspt = -*crosspt ; // flag this for second child

      pa = p1 ;
      pb = p2 ;
      } // If first child

   else {                       // Second child
      pa = p2 ;                 // so parents reverse roles
      pb = p1 ;
      } // If second child

/*
   Prepare for reproduction
*/

   if (*split) {              // Create left and right splitting masks
      right = 1 ;
      i = *split ;
      while (--i)
         right = (right << 1) | 1 ;
      left = 255 ^ right ;
      }

   if (*crosspt > 0) {        // Use one point crossover
      n1 = chromsize / 2 ;    // This many genes in first half of child
      n2 = chromsize - n1 ;   // and this many in second half
      n3 = n4 = 0 ;           // We are using one point crossover
      i = *crosspt - 1 ;      // We will start building child here
      }
   else {                             // Use two point crossover
      n1 = n2 = n3 = chromsize / 4 ;  // This many in first three quarters
      n4 = chromsize - n1 - n2 - n3 ; // And the last quarter gets the rest
      i = -*crosspt - 1 ;             // 2 point method was flagged by neg
      }

/*
   Do reproduction here
*/

   if (*split) {
      i = (i+1) % chromsize ;
      child[i] = (left & pa[i])  |  (right & pb[i]) ;
      --n1 ;
      }

   while (n1--) {
      i = (i+1) % chromsize ;
      child[i] = pb[i] ;
      }

   if (*split) {
      i = (i+1) % chromsize ;
      child[i] = (left & pb[i])  |  (right & pa[i]) ;
      --n2 ;
      }

   while (n2--) {
      i = (i+1) % chromsize ;
      child[i] = pa[i] ;
      }

   if (n4) {               // Two point crossover?

      if (*split) {
         i = (i+1) % chromsize ;
         child[i] = (left & pa[i])  |  (right & pb[i]) ;
         --n3 ;
         }

      while (n3--) {
         i = (i+1) % chromsize ;
         child[i] = pb[i] ;
         }

      if (*split) {
         i = (i+1) % chromsize ;
         child[i] = (left & pb[i])  |  (right & pa[i]) ;
         --n4 ;
         }

      while (n4--) {
         i = (i+1) % chromsize ;
         child[i] = pa[i] ;
         }

      } // If two point crossover
}

/*
--------------------------------------------------------------------------------

   mutate - apply the mutation operator to a single child

--------------------------------------------------------------------------------
*/
static void mutate (
   char *child ,   // Input/Output of the child
   int chromsize , // Number of variables in objective function
   double pmutate  // Probability of mutation
   )
{
   while (chromsize--) {
      if (unifrand() < pmutate)                          // Mutate this gene?
         child[chromsize] ^= (char) 1 << (longrand() % 8) ;  // Flip random bit
      }

}

/*
--------------------------------------------------------------------------------

   rand_ind - Randomly generate an individual's chromosome

--------------------------------------------------------------------------------
*/
static void rand_ind ( char *popptr , int chromsize )
{
   while (chromsize--)
      *popptr++ = 255 & longrand() ;
}

/*
--------------------------------------------------------------------------------

   decode - Decode the genes in this chromosome to the network weights.
            In genetic parlance, convert the genotype to the phenotype.

--------------------------------------------------------------------------------
*/

static void decode ( char *popptr , int nin , int nh1 , int nh2 ,
                     double *w1 , double *w2 )
{
   int n ;
   unsigned char gray, bit, parity, sum ;
   double *wptr ;

   if (! gray_initialized) {     // If the translation table has not yet
      gray_initialized = 1 ;     // been built, do it now (but just once!)
      for (n=0 ; n<256 ; n++) {  // Each gene is a one byte gray code
         gray = (unsigned char) n ;
         sum =  0 ;
         bit = 128 ;
         parity = 0 ;
         while (bit) {
            if (bit & gray)
               parity = ! parity ;
            if (parity)
               sum |= bit ;
            bit = bit >> 1 ;
            }
         gray_code_table[n] = sum ;
         }
      }

   n = nh1 * (nin+1) ; // Do hid layer 1 first.  It has this many weights.
   wptr = w1 ;         // Point to its weights

   while (n--) {
      gray = (unsigned char) *popptr++ ;  // The gene is a one byte gray code
      *wptr++ = ((double) gray_code_table[gray] - 127.5) * .0392 ; // -5 to 5
      }

   if (nh2) {              // Do second hidden layer if any
      n = nh2 * (nh1+1) ;
      wptr = w2 ;
      while (n--) {
         gray = (unsigned char) *popptr++ ;
         *wptr++ = ((double) gray_code_table[gray] - 127.5) * .0392 ;
         }
      }
}