de4_0.cpp

differential evolution algorithm optimization

/***************************************************************
**                                                            **
**            DE-Engine v4.0, Rainer Storn, 2004              **
**                                                            **
***************************************************************/

//------Compiler switches--------------------------------------------
#define DO_PLOTTING     //If defined, plotting is activated, and the functions in
                        //graphics.h are used (and might have to be adapted). Note
                        //that the entire project has to be compiled as a MS-Windows
                        //application.
                        //If not defined this is a simple console application and the project
                        //has to be compiled as such.   
//#define BOUND_CONSTR    //If defined the bounds fa_minbound[] and fa_maxbound[]
                        //are not only used for initializing the vector population
                        //but also used to keep the population within these bounds  

//------Include section----------------------------------------------

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
//#include <conio.h>    //not all C-Compilers have this (needed for function kbhit() )

#ifdef DO_PLOTTING
   #include "graphics.h"
#endif//DO_PLOTTING

#include "random.h"
#include "de.h"

//------Global variables--------------------------------------------

extern float gfa_bound[21];
#ifdef DO_PLOTTING
   float gfa_xmold[MAXDIM];  // "old" parameter vector (needed for erasing the old graph)
#endif//DO_PLOTTING
int   gi_gen;             // generation counter
int   gi_strategy;        // chooses DE-strategy
long  gl_nfeval;          // number of function evaluations     
int   gi_D;               // Dimension of parameter vector
int   gi_NP;              // Number of population members
int   gi_genmax;          // Maximum number of generations
t_pop gta_pop[2*MAXPOP];  // the two populations are put into one array side by side.      
t_pop gt_best;            // current best population member
t_pop *gpta_old, *gpta_new, *gpta_swap;



//------General functions-----------------------------------------

t_pop  evaluate(int gi_D, t_pop t_tmp, long *gl_nfeval, t_pop *tpa_array, int i_NP);
int    left_vector_wins(t_pop trial, t_pop target);
void   devol(FILE *Fp_in, FILE *Fp_out);
void   sort (t_pop ary[], int len);
void   assigna2b(int D, float a[], float b[]);
void   permute(int ia_urn2[], int i_urn2_depth, int i_NP, int i_avoid);


/**C*F****************************************************************
**                                                                  
** Function       :void  MainThread () / void main()                                        
**                                                                  
** Author         :Rainer Storn                                     
**                                                                  
** Description    :main program. Its makeup depends on whether the 
**                 project is set out to be a Windows application
**                 or a simple console application.                 
**                                                                  
** Functions      :-                                                
**                                                                  
** Globals        :gfa_bound[21]    (I)   problem dependent data vector
**                 gi_gen           (O)   generation counter 
**                 gi_strategy      (I)   DE-strategy
**                 gl_nfeval        (O)   number of function evaluations
**                 gi_D             (I)   number of parameters
**                 gt_best          (O)   best population member
**                                               
** Parameters     :-  
**                                                                  
** Preconditions  :-                     
**                                                                  
** Postconditions :- 
**
** Return Value   :-                                            
**                                                                  
***C*F*E*************************************************************/
#ifdef DO_PLOTTING //this is a MS-Windows application
 void MainThread ()
 {

     FILE *Fp_in;
     FILE *Fp_out;

	 Fp_in   = fopen("in.dat","r");
     Fp_out  = fopen("out.dat","w");             

     //-----Initialization of graphics----------------------
     graphics_init();

     //---optimization--------------------------------------
	 devol(Fp_in,Fp_out);

	 fclose(Fp_in);
	 fclose(Fp_out);

	 //-----enable zooming in and out even when optimization is finished----
	 while(1)
	 {
        update_graphics(gt_best.fa_vector, gi_D, gfa_bound, gl_nfeval, gi_gen, gt_best.fa_cost[0],gi_strategy,gi_genmax);
		Sleep(50); // provide some time (50ms) for the eye to see the graphics
        if (gi_exit_flag == 1)
		{
           gi_exit_flag = 0;
		}
	 }
 }

#else //just do console output
 void main()
 {

     FILE *Fp_in;
     FILE *Fp_out;

	 Fp_in   = fopen("in.dat","r");
     Fp_out  = fopen("out.dat","w");             

     //---optimization--------------------------------------
	 devol(Fp_in,Fp_out);

	 fclose(Fp_in);
	 fclose(Fp_out);
 }

#endif//DO_PLOTTING

void  assigna2b(int i_D, float fa_a[], float fa_b[])
/**C*F****************************************************************
**                                                                  
** Function       :void  assigna2b(int i_D, float fa_a[], float fa_b[])                                        
**                                                                  
** Author         :Rainer Storn                                     
**                                                                  
** Description    :Assigns i_D-dimensional vector fa_a to vector f_b.                 
**                                                                  
** Functions      :-                                                
**                                                                  
** Globals        :-
**                                               
** Parameters     :i_D     (I)     size of vectors
**                 fa_a[]  (I)     source vector
**                 fa_b[]  (I)     destination vector   
**                                                                  
** Preconditions  :-                     
**                                                                  
** Postconditions :- 
**
** Return Value   :-                                            
**                                                                  
***C*F*E*************************************************************/
{
   int j;
   for (j=0; j<i_D; j++)
   {
      fa_b[j] = fa_a[j];
   }
}


/**C*F****************************************************************
**                                                                  
** Function       :void devol(FILE *Fp_in, FILE *Fp_out)                                        
**                                                                  
** Author         :Rainer Storn, Kenneth Price                                     
**                                                                  
** Description    :Performs Differential Evolution optimization.                 
**                                                                  
** Functions      :-                                                
**                                                                  
** Globals        :
**                 gi_strategy  (I) 1 --> DE/rand/1:
**                                        the classical version of DE.             
**                                  2 --> DE/local-to-best/1:
**                                        a version which has been used by quite a number
**                                        of scientists. Attempts a balance between robustness
**                                        and fast convergence.
**                                  3 --> DE/best/1 with jitter:
**                                        taylored for small population sizes and fast convergence.
**                                        Dimensionality should not be too high.           
**                                  4 --> DE/rand/1 with per-vector-dither:
**                                        classical DE with dither to become even more robust.
**                                  5 --> DE/rand/1 with per-generation-dither:
**                                        classical DE with dither to become even more robust.
**                                        Choosing f_weight = 0.3 is a good start here.
**                                  6 --> DE/rand/1 either-or-algorithm:
**                                        Alternates between differential mutation and three-point-
**                                        recombination.           
**
**                 gi_D         (I)    number of parameters 
**                 gl_nfeval   (I/O)   counter for function evaluations                                               
**                 gi_gen      (I/O)   number of generations  
**                 gt_best     (I/O)   best vector
**                                               
** Parameters     :Fp_in    (I)     pointer to input file
**                 Fp_out   (I)     pointer to output file   
**                                                                  
** Preconditions  :-                     
**                                                                  
** Postconditions :- 
**
** Return Value   :-                                            
**                                                                  
***C*F*E*************************************************************/
void devol(FILE *Fp_in, FILE *Fp_out)
{
#define URN_DEPTH   5   //4 + one index to avoid

//-----Variable declarations-----------------------------------------

   int   i, j, k;                 // counting variables                 
   int   i_r1, i_r2, i_r3, i_r4;  // placeholders for random indexes    
   int   i_refresh;               // refresh rate of screen output      
   int   i_genmax, i_seed, i_bs_flag;
   int   ia_urn2[URN_DEPTH];

   float fa_minbound[MAXDIM];
   float fa_maxbound[MAXDIM];
   t_pop t_tmp, t_bestit;
#ifdef BOUND_CONSTR
   t_pop t_origin;
#endif//BOUND_CONSTR
   float f_weight, f_jitter, f_dither;
   float f_cross;


//-----Initialization of annealing parameters-------------------------

 fscanf(Fp_in,"%d",&gi_strategy);      //---choice of strategy-----------------
 fscanf(Fp_in,"%d",&i_genmax);         //---maximum number of generations------
 fscanf(Fp_in,"%d",&i_refresh);        //---output refresh cycle---------------
 gi_genmax = i_genmax;
 fscanf(Fp_in,"%d",&gi_D);             //---number of parameters---------------
 if (gi_D > MAXDIM)
 {
	 printf("Error! too many parameters\n");
	 return;
 }
 fscanf(Fp_in,"%d",&gi_NP);             //---population size.-------------------
 if (gi_NP > MAXPOP)
 {
	 printf("Error! too many points\n");
	 return;
 }
 fscanf(Fp_in,"%f",&f_weight);         //---weight factor----------------------
 fscanf(Fp_in,"%f",&f_cross);          //---crossing over factor---------------
 fscanf(Fp_in,"%d",&i_seed);           //---random seed------------------------
 fscanf(Fp_in,"%d",&i_bs_flag);        //---if TRUE: best of parent+child selection--------
                                       //---if FALSE: DE standard tournament selection-----
 for (i=0; i<gi_D; i++)
 {
	 fscanf(Fp_in,"%f",&fa_minbound[i]);
 }
 for (i=0; i<gi_D; i++)
 {
	 fscanf(Fp_in,"%f",&fa_maxbound[i]);
 }


//-----Initialize random number generator-----------------------------

 sgenrand((unsigned long)i_seed);// for Mersenne Twister
 gl_nfeval    =  0;              // reset number of function evaluations 

//------Initialization-----------------------------
   for (j=0; j<gi_D; j++)
   {
	 gta_pop[0].fa_vector[j] = fa_minbound[j]+genrand()*(fa_maxbound[j] - fa_minbound[j]);
   }
   gta_pop[0]      = evaluate(gi_D,gta_pop[0],&gl_nfeval,&gta_pop[0],gi_NP);
   gt_best  = gta_pop[0];

   for (i=1; i<gi_NP; i++)
   {
      for (j=0; j<gi_D; j++)
      {
	     gta_pop[i].fa_vector[j] = fa_minbound[j]+genrand()*(fa_maxbound[j] - fa_minbound[j]);
      }
      gta_pop[i] = evaluate(gi_D,gta_pop[i],&gl_nfeval,&gta_pop[0],gi_NP);

	  if (left_vector_wins(gta_pop[i],gt_best) == TRUE)
	  {
         gt_best = gta_pop[i];
	  }
   }

   t_bestit  = gt_best;

   //---assign pointers to current ("old") and new population---

   gpta_old = &gta_pop[0];
   gpta_new = &gta_pop[gi_NP];


//------Iteration loop--------------------------------------------
   gi_gen = 0;
#ifdef DO_PLOTTING
   while ((gi_gen < i_genmax) && (gi_exit_flag == 0))// && (gt_best.fa_cost[0] > VTR))