de4_0.cpp

differential evolution algorithm optimization

#else
   //Note that kbhit() needs conio.h which is not always available under Unix.
   while ((gi_gen < i_genmax))// && (kbhit() == 0))// && (gt_best.fa_cost[0] > VTR))
#endif//DO_PLOTTING
   {
      gi_gen++;

	  //----computer dithering factor (if needed)-----------------
	  f_dither = f_weight + genrand()*(1.0 - f_weight);

      //----start of loop through ensemble------------------------
      for (i=0; i<gi_NP; i++)           
      {
		permute(ia_urn2,URN_DEPTH,gi_NP,i); //Pick 4 random and distinct
		i_r1 = ia_urn2[1];                 //population members
		i_r2 = ia_urn2[2];
		i_r3 = ia_urn2[3];
		i_r4 = ia_urn2[4];

/*
        //---this is an alternative way to pick population members---
	    do                        // Pick a random population member 
		{
	       i_r1 = (int)(genrand()*gi_NP);
		}while(i_r1==i);

	    do                        // Pick a random population member 
		{
	       i_r2 = (int)(genrand()*gi_NP);
		}while((i_r2==i) || (i_r2==i_r1));

	    do                        // Pick a random population member 
		{
	       i_r3 = (int)(genrand()*gi_NP);
		}while((i_r3==i) || (i_r3==i_r1) || (i_r3==i_r2));

	    do                        // Pick a random population member 
		{
	       i_r4 = (int)(genrand()*gi_NP);
		}while((i_r4==i) || (i_r4==i_r1) || (i_r4==i_r2) || (i_r4==i_r3));
*/

        //========Choice of strategy=======================================================
		//---classical strategy DE/rand/1/bin-----------------------------------------
		if (gi_strategy == 1)
		{
           assigna2b(gi_D,gpta_old[i].fa_vector,t_tmp.fa_vector);
	       j = (int)(genrand()*gi_D); // random parameter         
	       k = 0;
	       do
		   {                            // add fluctuation to random target 
	          t_tmp.fa_vector[j] = gpta_old[i_r1].fa_vector[j] + f_weight*(gpta_old[i_r2].fa_vector[j]-gpta_old[i_r3].fa_vector[j]);

	          j = (j+1)%gi_D;
	          k++;
		   }while((genrand() < f_cross) && (k < gi_D));
#ifdef BOUND_CONSTR
		   assigna2b(gi_D,gpta_old[i_r1].fa_vector,t_origin.fa_vector);
#endif//BOUND_CONSTR
		}
		//---DE/local-to-best/1/bin---------------------------------------------------
	    else if (gi_strategy == 2)
		{
           assigna2b(gi_D,gpta_old[i].fa_vector,t_tmp.fa_vector);
	       j = (int)(genrand()*gi_D); // random parameter         
	       k = 0;
	       do
		   {                            // add fluctuation to random target 
	          t_tmp.fa_vector[j] = t_tmp.fa_vector[j] + f_weight*(t_bestit.fa_vector[j] - t_tmp.fa_vector[j]) + 
				                                        f_weight*(gpta_old[i_r2].fa_vector[j]-gpta_old[i_r3].fa_vector[j]);

	          j = (j+1)%gi_D;
	          k++;
		   }while((genrand() < f_cross) && (k < gi_D));
#ifdef BOUND_CONSTR
		   assigna2b(gi_D,t_tmp.fa_vector,t_origin.fa_vector);
#endif//BOUND_CONSTR
		}
		//---DE/best/1/bin with jitter------------------------------------------------
	    else if (gi_strategy == 3)
		{
           assigna2b(gi_D,gpta_old[i].fa_vector,t_tmp.fa_vector);
	       j = (int)(genrand()*gi_D); // random parameter         
	       k = 0;
	       do
		   {                            // add fluctuation to random target 
			  f_jitter = (0.0001*genrand()+f_weight);
	          t_tmp.fa_vector[j] = t_bestit.fa_vector[j] + f_jitter*(gpta_old[i_r1].fa_vector[j]-gpta_old[i_r2].fa_vector[j]);

	          j = (j+1)%gi_D;
	          k++;
		   }while((genrand() < f_cross) && (k < gi_D));
#ifdef BOUND_CONSTR
		   assigna2b(gi_D,t_tmp.fa_vector,t_origin.fa_vector);
#endif//BOUND_CONSTR
		}
		//---DE/rand/1/bin with per-vector-dither-------------------------------------
	    else if (gi_strategy == 4)
		{
           assigna2b(gi_D,gpta_old[i].fa_vector,t_tmp.fa_vector);
	       j = (int)(genrand()*gi_D); // random parameter         
	       k = 0;
	       do
		   {                            // add fluctuation to random target 
	          t_tmp.fa_vector[j] = gpta_old[i_r1].fa_vector[j] + 
			                       (f_weight + genrand()*(1.0 - f_weight))*
								   (gpta_old[i_r2].fa_vector[j]-gpta_old[i_r3].fa_vector[j]);

	          j = (j+1)%gi_D;
	          k++;
		   }while((genrand() < f_cross) && (k < gi_D));
#ifdef BOUND_CONSTR
		   assigna2b(gi_D,t_tmp.fa_vector,t_origin.fa_vector);
#endif//BOUND_CONSTR
		}
		//---DE/rand/1/bin with per-generation-dither---------------------------------
	    else if (gi_strategy == 5)
		{
           assigna2b(gi_D,gpta_old[i].fa_vector,t_tmp.fa_vector);
	       j = (int)(genrand()*gi_D); // random parameter         
	       k = 0;
	       do
		   {                            // add fluctuation to random target 
	          t_tmp.fa_vector[j] = gpta_old[i_r1].fa_vector[j] + f_dither*(gpta_old[i_r2].fa_vector[j]-gpta_old[i_r3].fa_vector[j]);

	          j = (j+1)%gi_D;
	          k++;
		   }while((genrand() < f_cross) && (k < gi_D));
#ifdef BOUND_CONSTR
		   assigna2b(gi_D,t_tmp.fa_vector,t_origin.fa_vector);
#endif//BOUND_CONSTR
		}
		//---variation to DE/rand/1/bin: either-or-algorithm--------------------------
	    else
		{
           assigna2b(gi_D,gpta_old[i].fa_vector,t_tmp.fa_vector);
	       j = (int)(genrand()*gi_D); // random parameter         
	       k = 0;
		   if (genrand() < 0.5) //Pmu = 0.5
		   {//differential mutation
	          do
			  {                            // add fluctuation to random target 
	             t_tmp.fa_vector[j] = gpta_old[i_r1].fa_vector[j] + f_weight*(gpta_old[i_r2].fa_vector[j]-gpta_old[i_r3].fa_vector[j]);

	             j = (j+1)%gi_D;
	             k++;
			  }while((genrand() < f_cross) && (k < gi_D));
		   }
		   else
		   {//recombination with K = 0.5*(F+1) --> F-K-Rule
	          do
			  {                            // add fluctuation to random target 
	             t_tmp.fa_vector[j] = gpta_old[i_r1].fa_vector[j] + 0.5*(f_weight+1.0)*
					                 (gpta_old[i_r2].fa_vector[j]+gpta_old[i_r3].fa_vector[j] -
									2*gpta_old[i_r1].fa_vector[j]);

	             j = (j+1)%gi_D;
	             k++;
			  }while((genrand() < f_cross) && (k < gi_D));
		   }
#ifdef BOUND_CONSTR
		   assigna2b(gi_D,gpta_old[i_r1].fa_vector,t_origin.fa_vector);
#endif//BOUND_CONSTR		
		}//end if (gi_strategy ...

#ifdef BOUND_CONSTR
      for (j=0; j<gi_D; j++) //----boundary constraints via random reinitialization-------
      {                      //----and bounce back----------------------------------------
         if (t_tmp.fa_vector[j] < fa_minbound[j])
		 {
            t_tmp.fa_vector[j] = fa_minbound[j]+genrand()*(t_origin.fa_vector[j] - fa_minbound[j]);
		 }
         if (t_tmp.fa_vector[j] > fa_maxbound[j])
		 {
            t_tmp.fa_vector[j] = fa_maxbound[j]+genrand()*(t_origin.fa_vector[j] - fa_maxbound[j]);
		 }
      }
#endif//BOUND_CONSTR


        //------Trial mutation now in t_tmp-----------------

	    t_tmp = evaluate(gi_D,t_tmp,&gl_nfeval,&gpta_old[0],gi_NP);  // Evaluate mutant in t_tmp[] 

if (i_bs_flag == TRUE)
{
        gpta_new[i]=t_tmp; //save new vector, selection will come later
}
else
{
		if (left_vector_wins(t_tmp,gpta_old[i]) == TRUE)
		{
	       gpta_new[i]=t_tmp;              // replace target with mutant 

		   if (left_vector_wins(t_tmp,gt_best) == TRUE)// Was this a new minimum? 
		   {                               // if so...
		      gt_best = t_tmp;             // store best member so far  
		   }                               // If mutant fails the test...
		}                                  // go to next the configuration 
	    else
		{
	       gpta_new[i]=gpta_old[i];              // replace target with old value 
		}
}//if (i_bs_flag == TRUE)
      }//end for (i=0; i<gi_NP; i++)
					    // End mutation loop through pop. 

if (i_bs_flag == TRUE)
{
      sort (gpta_old, 2*gi_NP); //sort array of parents + children
	  gt_best = gpta_old[0];
}
else
{
      gpta_swap = gpta_old;
      gpta_old  = gpta_new;
      gpta_new  = gpta_swap;
}//if (i_bs_flag == TRUE)

	  t_bestit = gt_best;


//    }//if ....
//======Output Part=====================================================

	if (gi_gen%i_refresh == 0) //refresh control
	{
       #ifdef DO_PLOTTING
          update_graphics(gt_best.fa_vector, gi_D, gfa_bound, gl_nfeval, gi_gen, gt_best.fa_cost[0], gi_strategy,gi_genmax);
       #else
	      printf("%6d   %12.6f   %12.6f\n",gl_nfeval, gt_best.fa_cost[0], gt_best.fa_constraint[0]);
       #endif//DO_PLOTTING
	   fprintf(Fp_out,"%6d   %12.6f\n",gl_nfeval, gt_best.fa_cost[0]);


	}//end if (gi_gen%10 == 0)

   }//end while ((gi_gen < i_genmax) && (gf_emin > MINI))

   fprintf(Fp_out,"\n******** best vector ********\n", i, gt_best.fa_vector[i]);
   for (i=0; i<gi_D; i++)
   {
	   #ifndef DO_PLOTTING
	      printf("%12.6f\n", gt_best.fa_vector[i]);
	   #endif//DO_PLOTTING
	   fprintf(Fp_out,"best_vector[%d]=%12.6f\n", i, gt_best.fa_vector[i]);
   }
   #ifdef DO_PLOTTING 
      if (gi_exit_flag == 1)
	  {
         gi_exit_flag = 0;
	  }
   #endif//DO_PLOTTING
}


void permute(int ia_urn2[], int i_urn2_depth, int i_NP, int i_avoid)
/**C*F****************************************************************
**                                                                  
** Function       :void permute(int ia_urn2[], int i_urn2_depth)                                        
**                                                                  
** Author         :Rainer Storn                                     
**                                                                  
** Description    :Generates i_urn2_depth random indices ex [0, i_NP-1]
**                 which are all distinct. This is done by using a 
**                 permutation algorithm called the "urn algorithm"
**                 which goes back to C.L.Robinson.                
**                                                                  
** Functions      :-                                                
**                                                                  
** Globals        :-
**                                               
** Parameters     :ia_urn2       (O)    array containing the random indices
**                 i_urn2_depth  (I)    number of random indices (avoided index included)
**                 i_NP          (I)    range of indices is [0, i_NP-1]
**                 i_avoid       (I)    is the index to avoid and is located in
**                                      ia_urn2[0].   
**                                                                  
** Preconditions  :# Make sure that ia_urn2[] has a length of i_urn2_depth.
**                 # i_urn2_depth must be smaller than i_NP.                     
**                                                                  
** Postconditions :# the index to be avoided is in ia_urn2[0], so fetch the
**                   indices from ia_urn2[i], i = 1, 2, 3, ..., i_urn2_depth. 
**
** Return Value   :-                                            
**                                                                  
***C*F*E*************************************************************/

{
	int  i, k, i_urn1, i_urn2;
	int  ia_urn1[MAXPOP] = {0};      //urn holding all indices

	k      = i_NP;
	i_urn1 = 0; 
	i_urn2 = 0;
	for (i=0; i<i_NP; i++) ia_urn1[i] = i; //initialize urn1

	i_urn1 = i_avoid;                  //get rid of the index to be avoided and place it in position 0.
	while (k >= i_NP-i_urn2_depth)     //i_urn2_depth is the amount of indices wanted (must be <= NP) 
	{
	   ia_urn2[i_urn2] = ia_urn1[i_urn1];      //move it into urn2
	   ia_urn1[i_urn1] = ia_urn1[k-1]; //move highest index to fill gap
	   k = k-1;                        //reduce number of accessible indices
	   i_urn2 = i_urn2 + 1;            //next position in urn2
       i_urn1 = (int)(genrand()*k);    //choose a random index
	}
}