main.c

用C编写的部落寻优的优化算法

					}
					i=i-1;
				}
				TR[level].size=TR[level].size-1;
				if (problem[level].printlevel>0)

				{
					printf("\n %i tribes. No particle anymore in tribe %i!",TR[level].size+1,i0);
					printf(" => New tribe number %i",TR[level].size);
				}
				if (TR[level].size==0)
				{
					printf("\n NO PARTICLE ANYMORE!");
					goto the_end;
				}
				goto empty_tribe;
		}

if(linkreorg==1 && TR[level].size>1) // NEW 6.2 Optionnally reorganize links
{
	link_reorg(level);
}

		cycle=0;
		if(problem[level].save==1) save_swarm(f_swarm,level);

	end_adapt:;       //   ---------------------------  ADAPTATION(S)  (end)



 if(rand_hood>0)
 {
      swarm_size=problem[level].init_size+n_add-n_remove;
	  if(circular_hood==1) rand_hood=K_optim(swarm_size,problem[level].eps);
      for(k=0;k<TR[level].size;k++)
      {
         size= TR[level].tr[k].size;
         for (i=0;i<size;i++)
         {
            TR[level].tr[k].part[i].ig.size=rand_hood;

            for (j=0;j<TR[level].tr[k].part[i].ig.size;j++)
            {
               rank=alea(0,swarm_size);
               TR[level].tr[k].part[i].ig.label[j]=TR[level].tr[k].part[rank].label;
            }
         }
      }
 }

         if (problem[level].save==-1)  // Save swarm size and energies
										// Just for information
        {
           Ek=energy_k(level); Ep=energy_p(level);
		   swarm_size=problem[level].init_size+n_add-n_remove;
           fprintf(f_energy,"\n%i %i %i %i %f %f",(int)eval_f_tot[level],swarm_size,n_add,n_remove,Ek,Ep);
        }


		goto do_Iter; // Iterates till a criterion is met

end_by_chance:
		printf("\n CHANCE !");

the_end: // ------------------------------------------------------------------------------------------------------------------- END
	error_cont=offline_error_cont/(eval_f_tot[level]-n_add);
	if (DYN==0)
    error=total_error(best_result.f);
	else error=offline_error/n_change;// Mean error for dynamic optimisation

	if (problem[level].save==-2) // Save info for convergence curve plotting
	{
		fprintf(f_trace_run,"\n%.0f %f %f %i",eval_f_tot[level],error,error_cont,n_add);
	}

ticks2=clock();
clock_tick=ticks2-ticks1; // Processor time needed
duration = (double)(clock_tick/CLOCKS_PER_SEC);
duration_tot=duration_tot+duration;

eval_run[times]=eval_f_tot[level];
eps_run[times][0]=error;
eps_run[times][1]=error_cont;

			// Keep the min. and max. eps. Just for information
			if (error<min_eps) min_eps=error;
			if (error>max_eps) max_eps=error;

if (level>0)  goto end_PSO;

mean_swarm_size=mean_swarm_size/Iter;
Mean_swarm=Mean_swarm+mean_swarm_size;

if (problem[level].printlevel>=0)
{
   if(error<=problem[level].eps)
   {
      printf("\n SUCCESS");
      }
      else printf("\n FAILURE");

  printf ("\n Run %i => eps %f, eps_cont %f,eval: %.0f, duration: %.2f",times+1,error,error_cont,eval_f_tot[level],duration);
  printf(" + %i - %i",n_add,n_remove);
  printf(" (%.0f)",mean_swarm_size);
   display_position(best_result);



//fprintf(f_trace,"%f %f\n",error,eval_f_tot[level]);


// Save run result
// This print should be consistent with the titles written at the very beginning
fprintf (f_run,"\n%i %i %f",problem[level].funct,times+1,problem[level].target);
if(DYN==0)
for (i=0;i<best_result.f.size;i++) fprintf (f_run," %f",best_result.f.f[i]); // f value(s)
else fprintf(f_run," %.12f",error);// For non recursive dynamic optimisation

fprintf(f_run," %.0f %f %i %i",eval_f_tot[level],duration,n_add,n_remove);

for (i=0;i<best_result.p.size;i++) fprintf(f_run," %.12f",best_result.p.x[i]); // position

  if (problem[level].funct==11) // Apple trees
  {
    post= homogen_to_carte(best_result);
    printf("\n   Cartesian coordinates: ");
    for (i=0;i<post.p.size;i++)  printf(" %f",post.p.x[i]);
      fprintf(f_run,"   ");
  for (i=0;i<post.p.size;i++)  fprintf(f_run," %f",post.p.x[i]);
  }
}


if (eval_f_tot[level]>nb_eval_max) nb_eval_max=eval_f_tot[level];
if (eval_f_tot[level]<nb_eval_min) nb_eval_min=eval_f_tot[level];


if(problem[level].printlevel>0) print_N_run(level,best_result);

if(error<total_error(BEST.f)) BEST=best_result; // Best of the best  amongst runs

	times=times+1;
	if (times<problem[level].times) goto times_loop;

  //-----------------------------------------------------------------------------------------
  if (problem[level].printlevel<0) goto end_PSO;

		ticks2=clock();
		mean_eval=zzz/problem[level].times;
		mean_eps=mean_eps/problem[level].times;
		Mean_swarm=Mean_swarm/problem[level].times;
		duration=duration_tot/problem[level].times;

		mean_eval=0;
		mean_eps=0; mean_eps_cont=0;
//		times=problem[level].times;
		for (i=0;i<times;i++)
		{
			mean_eval=mean_eval+eval_run[i];
			mean_eps=mean_eps+eps_run[i][0];
			mean_eps_cont=mean_eps_cont+eps_run[i][1];
		}
		mean_eval=mean_eval/times;
		mean_eps=mean_eps/times;
		mean_eps_cont=mean_eps_cont/times;

		sigma_eval=0;
		for (i=0;i<problem[level].times;i++)
		{
			sigma_eval=sigma_eval+(eval_run[i]-mean_eval)*(eval_run[i]-mean_eval);
		}
		sigma_eval=sqrt (sigma_eval/problem[level].times);


		sigma_eps=0; sigma_eps_cont=0;
		for (i=0;i<problem[level].times;i++)

		{
			sigma_eps=sigma_eps+(eps_run[i][0]-mean_eps)*(eps_run[i][0]-mean_eps);
			sigma_eps_cont+=(eps_run[i][1]-mean_eps_cont)*(eps_run[i][1]-mean_eps_cont);
		}
		sigma_eps=sqrt (sigma_eps/problem[level].times);
		sigma_eps_cont=sqrt (sigma_eps_cont/problem[level].times);

		fprintf(f_synth,"\n");
		fprintf(f_synth,"%2i",problem[level].funct);
		fprintf(f_synth," %2i",problem[level].DD);
		fprintf(f_synth," %4.3f",problem[level].target);
		fprintf(f_synth," %4.10f",problem[level].eps);
		fprintf(f_synth," %i",problem[level].times);
		fprintf(f_synth," %.0f",mean_eval); // Mean number of evaluations
		fprintf(f_synth," %.0f",sigma_eval); // Standard deviation for eval
      z= 1- (double)failure_tot/(double)problem[level].times ;
		fprintf(f_synth," %.3f",z);  // Success rate

		fprintf(f_synth, " %.0f %.0f",nb_eval_min,nb_eval_max);
		fprintf(f_synth, " %4.0f ",Mean_swarm);

		fprintf(f_synth, " %4.10f ",mean_eps);
		fprintf(f_synth," %4.10f",sigma_eps); // Standard deviation for eps
		fprintf(f_synth, " %4.10f ",mean_eps_cont);
		fprintf(f_synth," %4.10f",sigma_eps_cont);

		fprintf(f_synth, " %4.10f ",min_eps);
		fprintf(f_synth, " %4.10f ",max_eps);

		fprintf(f_synth," / Time %4.4f",duration);

		printf("\nMEAN NUMBER OF EVALUATIONS: %.0f",mean_eval); // Mean number of evaluations
		printf(" sigma %.0f",sigma_eval); // Standard deviation
		printf(" [%.0f %.0f ]",nb_eval_min,nb_eval_max);

		printf("\nMean eps: %4.10f",mean_eps);
		printf(" sigma %4.10f",sigma_eps);
		printf(" [%4.10f, %4.10f]",min_eps,max_eps);
		printf("\nMean eps cont.: %4.10f",mean_eps_cont);
		printf(" sigma %4.10f",sigma_eps_cont);


		printf("\nMean swarm:%4.0f ",Mean_swarm);

		printf("\nTOTAL NUMBER OF FAILURES: %i/%i", failure_tot,problem[level].times);
//		if (eval_f_sup>0) printf(" (after %.0f evaluations)",eval_f_sup);
  printf("\n Success rate %.3f", z );


		duration = duration_tot/(float)problem[level].times;
		printf(" / MEAN TIME %4.4f  ",duration);

/*
   if (level==0)    // Display strategies used
   {
      printf("\n strategies used: ");
      for (i=0;i<nb_strateg;i++)
      if (strateg[i]>0)

         printf("\n %i: %.0f times  ",i,strateg[i]) ;
     }
 */

 if (level==0) // Display particle status
 {
   printf("\n\n Status  Count    Rate ");
    z=0; for (i=0;i<Max_status;i++) z=z+status_count[i];
   for (i=0;i<Max_status;i++) printf("\n %i      %.0f    %.0f",i,status_count[i],100.0*status_count[i]/z);
 }

 if (problem[level].printlevel>1) display_swarm(1,level);

 save_swarm(f_swarm,level);
fprintf(f_swarm,"\n%i\n",-1); // End of save swarm file
      fclose(f_swarm);
         fprintf(f_run,"\n%i\n",-1); // End of save run file
         fclose(f_run);

  if (problem[level].nb_f==1) goto end_PSO;

  // Multiobjective problem. Clean up the result file
printf("\n Multiobjective cleaning");
  f_run=fopen("run.txt","r");
  f_run_clean=fopen("run_cleaned.txt","w");

  clean_run(f_run,f_run_clean,problem[level].nb_f,problem[level].DD,level);

end_PSO:
if (problem[level].printlevel>0)
{
   if(MEMO==1) display_memo(level);
}
return best_result;

}


 // ----------------------------------------------------------------------------- ADD_MEMO
void add_memo(struct position x, int level)
{
 /*
   Put a new position in memo, so that memo is sorted by increasing error.
   Note: memo[2] is a global variable
   Useful only if MEMO option is used.
 */

 double error,error1;
  int  i;
  int rank_p;

 error=total_error(x.f);
 if (memo[level].size==0)  // The very first time
 {
       memo[level].x[memo[level].size]=x;
       memo[level].error[memo[level].size]=error;
       memo[level].size=memo[level].size+1;
       rank_p=0;
       goto end;
 }


 // Find the first worst already stored
 for (i=0;i<memo[level].size;i++)
 {
   error1=memo[level].error[i];
   if (error1>=error) {rank_p=i;goto shift;}
 }

 rank_p=memo[level].size;
 goto insert;

 shift:  // Make place for the new position to store

  for (i=memo[level].size;i>rank_p;i=i-1)
  {
     memo[level].x[i]=  memo[level].x[i-1] ;
     memo[level].error[i]=memo[level].error[i-1];
  }



 insert:  // Put the new position
   memo[level].size=memo[level].size+1;
  if (memo[level].size>=Max_Memo-1) memo[level].size=Max_Memo-1;

   memo[level].x[rank_p]=x;
   memo[level].error[rank_p]=error;

   end:
 //printf("\n add_memo (size %i) rank %i, error %f",memo[level].size,rank_p,error);
 //display_memo(level)
;

}

// ----------------------------------------------------------------------------- BEST_INFORMER
struct particle	best_informer(struct particle part,int no_best,int level)
{
/*
   Look for the best informer of the particle.
   If no_best=1, choose it at random.
    2 => pseudo gradient
	0 (default) => really the best
	3 => 0 or 2 at random


*/
int               better;
char			bid;
double			delta1,delta2;
double			dist;
double			error;
double			grad;
int				i;
int				j;

int				j_best;
int				k;
int				k_best;

struct f			f0;
int				labelt;

switch (no_best)
{
case 1:// Choose it at random

	i=alea(0,part.ig.size-1);
    labelt=part.ig.label[i];
	// Find the particle with this label
	for (j=0;j<TR[level].size;j++) // For each tribe ...
		for (k=0;k<TR[level].tr[j].size;k++) // ... for each particle in the tribe
		{
			if (TR[level].tr[j].part[k].label==labelt)
               {
                j_best=j;
                k_best=k;
                break;
               }
		}
	printf("\nERROR. best_informer 1. Can't find particle with label %i",labelt);
	display_swarm(1,level);
	scanf("%c",&bid);
    return TR[level].tr[j].part[k];

case 2: // Pseudo gradient
case2:
	j_best=-1;
	k_best=0;
	grad=0;

	for (i=0;i<part.ig.size;i++) // For each informer
	{
		labelt=part.ig.label[i];

		// Find the particle with this label
		for (j=0;j<TR[level].size;j++) // For each tribe ...
		{
			for (k=0;k<TR[level].tr[j].size;k++) // ... for each particle in the tribe
			{
				if (TR[level].tr[j].part[k].label==labelt) goto compare2;
			}
		}
		printf("\nERROR. best_informer (2). Can't find particle with label %i",labelt);
		printf("\n Particle label %i. ig list:",part.label);
		for (k=0;k<part.ig.size;k++) printf(" %i",part.ig.label[k]);
		display_swarm(1,level);
		scanf("%c",&bid);

		break;


	compare2:
		error=total_error(part.p_i.f);
		delta1=error-total_error(TR[level].tr[j].part[k].p_i.f);
		dist=distance(part.p_i,TR[level].tr[j].part[k].p_i);
		if (dist<=0) continue;


		delta2=delta1/dist;
		if (delta2>grad)
		{
			j_best=j;k_best=k; grad=delta2;
		}
	}
  if(j_best>=0) break;
  // If no better than the particle itself
   //if(j_best<0) return part;  // return the part
   if(j_best<0) goto case0; // Use the classical method
   if(j_best<0) {j_best=0;k_best=0;break;} // Choose particle 0


case 3: // Random 0 or 2

	if (alea(0,1)<0.5) goto case0; else goto case2;

default: // Find really the best
case0:
	j_best=0;
	k_best=0;
	//f0=TR[level].tr[j_best].part[k_best].p_i.f;
	f0.size=TR[level].tr[j_best].part[k_best].p_i.f.size;
	for (i=0;i<f0.size;i++) f0.f[i]=infinite;