ga_simplex.c

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

 * Evaluate the function at this reflected point.  
 */
    pop->simplex_params->from_double(pop, new1, new1_d);
    pop->evaluate(pop, new1);

    if (new1->fitness > putative[0]->fitness)
      {
/*
 * The new solution is fitter than the previously fittest solution, so attempt an 
 * additional extrapolation by a factor alpha.
 */
/*printf("DEBUG: new1 (%f) is fitter than p0 ( %f )\n", new1->fitness, putative[0]->fitness);*/

      for (j = 0; j < pop->simplex_params->dimensions; j++)
        new2_d[j] = (1.0 + pop->simplex_params->alpha) * new1_d[j] -
                    pop->simplex_params->alpha * putative_d[num_points-1][j];

      pop->simplex_params->from_double(pop, new2, new2_d);
      pop->evaluate(pop, new2);

      if (new2->fitness > putative[0]->fitness)
        {
/*
 * This additional extrapolation succeeded, so replace the least fit solution
 * by inserting new solution in correct position.
 */
/*printf("DEBUG: new2 (%f) is fitter than p0 ( %f )\n", new2->fitness, putative[0]->fitness);*/

        tmpentity = putative[pop->simplex_params->dimensions];
        tmpdoubleptr = putative_d[pop->simplex_params->dimensions];

        for (j = pop->simplex_params->dimensions; j > 0; j--)
          {
          putative[j]=putative[j-1];
          putative_d[j]=putative_d[j-1];
          }

        putative[0] = new2;
        putative_d[0] = new2_d;

        new2 = tmpentity;
        new2_d = tmpdoubleptr;
        }
      else
        {
/*
 * This additional extrapolation failed, so use the original
 * reflected solution.
 */
        tmpentity = putative[pop->simplex_params->dimensions];
        tmpdoubleptr = putative_d[pop->simplex_params->dimensions];

        for (j = pop->simplex_params->dimensions; j > 0; j--)
          {
          putative[j]=putative[j-1];
          putative_d[j]=putative_d[j-1];
          }

        putative[0] = new1;
        putative_d[0] = new1_d;

        new1 = tmpentity;
        new1_d = tmpdoubleptr;
        }
      }
    else if (new1->fitness < putative[pop->simplex_params->dimensions-1]->fitness)
      {
/*
 * The reflected point is worse than the second-least fit.  
 */
/*printf("DEBUG: new1 (%f) is less fit than p(n-1) ( %f )\n", new1->fitness, putative[pop->simplex_params->dimensions-1]->fitness);*/

      did_replace = FALSE;

      if (new1->fitness > putative[pop->simplex_params->dimensions]->fitness)
        {
/*
 * It is better than the least fit, so use it to replace the
 * least fit.
 */
/*printf("DEBUG: but fitter than p(n) ( %f )\n", putative[pop->simplex_params->dimensions]->fitness);*/
        did_replace = TRUE;

        tmpentity = putative[pop->simplex_params->dimensions];
        tmpdoubleptr = putative_d[pop->simplex_params->dimensions];

        putative[pop->simplex_params->dimensions] = new1;
        putative_d[pop->simplex_params->dimensions] = new1_d;

        new1 = tmpentity;
        new1_d = tmpdoubleptr;
        }
/*
 * Perform a contraction of the simplex along one dimension, away from worst point.
 */
      for (j = 0; j < pop->simplex_params->dimensions; j++)
        new1_d[j] = (1.0 - pop->simplex_params->beta) * average[j] +
                   pop->simplex_params->beta * putative_d[num_points-1][j];

      pop->simplex_params->from_double(pop, new1, new1_d);
      pop->evaluate(pop, new1);

      if (new1->fitness > putative[pop->simplex_params->dimensions]->fitness)
        {
/*
 * The contraction gave an improvement, so accept it by
 * inserting the new solution at the correct position.
 */
        did_replace = TRUE;

/*printf("DEBUG: contracted new1 (%f) is fitter than p(n) ( %f )\n", new1->fitness, putative[pop->simplex_params->dimensions]->fitness);*/

        i = 0;
        while (putative[i]->fitness > new1->fitness) i++;

        tmpentity = putative[pop->simplex_params->dimensions];
        tmpdoubleptr = putative_d[pop->simplex_params->dimensions];

        for (j = pop->simplex_params->dimensions; j > i; j--)
          {
          putative[j]=putative[j-1];
          putative_d[j]=putative_d[j-1];
          }

        putative[i] = new1;
        putative_d[i] = new1_d;

        new1 = tmpentity;
        new1_d = tmpdoubleptr;
        }

      if (did_replace == FALSE)
        {
/*
 * The new solution is worse than the previous worse.  So, contract
 * toward the average point.
 */
/*printf("DEBUG: new1 (%f) is worse than all.\n", new1->fitness);*/

        for (i = 1; i < num_points; i++)
          {
          for (j = 0; j < pop->simplex_params->dimensions; j++)
            putative_d[i][j] = average[j] +
                               pop->simplex_params->gamma * (putative_d[i][j] - average[j]);

          pop->simplex_params->from_double(pop, putative[i], putative_d[i]);
          pop->evaluate(pop, putative[i]);
          }

/*
 * Alternative is to contact toward the most fit point.
        for (i = 1; i < num_points; i++)
          {
          for (j = 0; j < pop->simplex_params->dimensions; j++)
            putative_d[i][j] = putative_d[0][j] +
                               pop->simplex_params->gamma * (putative_d[i][j] - putative_d[0][j]);

          pop->simplex_params->from_double(pop, putative[i], putative_d[i]);
          pop->evaluate(pop, putative[i]);
          }
*/

        }

      }
    else
      {
/*
 * The reflection gave a solution which was better than the worst two
 * solutions, but worse than the best solution.
 * Replace the old worst solution by inserting the new solution at the
 * correct position.
 */
/*printf("DEBUG: new1 (%f) is fitter than worst 2\n", new1->fitness);
      for (j=0; j < pop->simplex_params->dimensions; j++)
        printf("%d fitness = %f\n", j, putative[j]->fitness);
*/

      i = 0;
      while (putative[i]->fitness > new1->fitness) i++;

/*printf("DEBUG: new1 inserted at position %d\n", i);*/

      tmpentity = putative[pop->simplex_params->dimensions];
      tmpdoubleptr = putative_d[pop->simplex_params->dimensions];

      for (j = pop->simplex_params->dimensions; j > i; j--)
        {
        putative[j]=putative[j-1];
        putative_d[j]=putative_d[j-1];
        }

      putative[i] = new1;
      putative_d[i] = new1_d;

      new1 = tmpentity;
      new1_d = tmpdoubleptr;
      }

/*
 * Use the iteration callback.
 */
    plog( LOG_VERBOSE,
          "After iteration %d, the current solution has fitness score of %f",
          iteration,
          putative[0]->fitness );

    }	/* Iteration loop. */

/*
 * Store best solution.
 */
  ga_entity_copy(pop, initial, putative[0]);

/*
 * Cleanup.
 */
  ga_entity_dereference(pop, new1);
  ga_entity_dereference(pop, new2);

  for (i=0; i<num_points; i++)
    {
    ga_entity_dereference(pop, putative[i]);
    }

  s_free(putative);
  s_free(putative_d);
  s_free(putative_d_buffer);

  return iteration;
  }


/**********************************************************************
  ga_simplex_double()
  synopsis:	Performs optimisation on the passed entity by using a
  		simplistic simplex-search.  The fitness evaluations
		are performed using the standard and evaluation
		callback mechanism.
		The passed entity will have its data overwritten.  The
		remainder of the population will be let untouched.
		Note that it is safe to pass a NULL initial structure,
		in which case a random starting structure will be
		generated, however the final solution will not be
		available to the caller in any obvious way.

		Only double chromosomes may be used in this optimised
		version of the algorithm.
  parameters:
  return:
  last updated:	13 Apr 2005
 **********************************************************************/

int ga_simplex_double(	population		*pop,
		entity			*initial,
		const int		max_iterations )
  {
  int		iteration=0;		/* Current iteration number. */
  int		i, j;			/* Index into putative solution array. */
  entity	**putative;		/* Current working solutions. */
  entity	*new1, *new2;		/* New putative solutions. */
  entity	*tmpentity;		/* Used to swap working solutions. */
  int		num_points;		/* Number of search points. */
  double	*average;		/* Vector average of solutions. */
  int           first=0, last;		/* Indices into solution arrays. */
  boolean       done=FALSE;		/* Whether the shuffle sort is complete. */
  boolean	did_replace;		/* Whether worst solution was replaced. */
  boolean	restart_needed;		/* Whether the search needs restarting. */

/*
 * Checks.
 */
  if (!pop) die("NULL pointer to population structure passed.");
  if (!pop->evaluate) die("Population's evaluation callback is undefined.");
  if (!pop->simplex_params) die("ga_population_set_simplex_params(), or similar, must be used prior to ga_simplex().");

/* 
 * Prepare working entities and double arrays.
 * The space for the average and new arrays are allocated simultaneously.
 */
  num_points = pop->len_chromosomes+1;
  putative = s_malloc(sizeof(entity *)*num_points);
  average = s_malloc(sizeof(double)*pop->len_chromosomes);

  for (i=1; i<num_points; i++)
    {
    putative[i] = ga_get_free_entity(pop);    /* The 'working' solutions. */
    }

  new1 = ga_get_free_entity(pop);
  new2 = ga_get_free_entity(pop);

/* Do we need to generate a random starting solution? */
  if (!initial)
    {
    plog(LOG_VERBOSE, "Will perform simplex search with random starting solution.");

    putative[0] = ga_get_free_entity(pop);
    ga_entity_seed(pop, putative[0]);
    initial = ga_get_free_entity(pop);
    }
  else
    {   
    plog(LOG_VERBOSE, "Will perform simplex search with specified starting solution.");

    putative[0] = ga_get_free_entity(pop);
    ga_entity_copy(pop, putative[0], initial);
    }

/*
 * Generate sample points.
 * Ensure that these initial solutions are scored.
 *
 * NOTE: Only perturb each solution by one dimension, by an
 * amount specified by the step parameter; it might be better to perturb
 * all dimensions and/or by a randomized amount.
 */
#pragma omp parallel \
   if (GAUL_DETERMINISTIC_OPENMP==0) \
   shared(pop,num_points,putative) private(i,j)
    {
#pragma omp single \
   nowait
    pop->evaluate(pop, putative[0]);

#pragma omp for \
   schedule(static) nowait
    for (i=1; i<num_points; i++)