pso8.c

C语言开发的微粒群优化算法源程序

		if (level==0)
		{
			fprintf(f_synth,"\n");
			fprintf(f_synth," %2i",param.funct);
			fprintf(f_synth," %2i",param.DD);
			fprintf(f_synth," %4.3f",param.target);
			fprintf(f_synth," %4.10f",param.eps);
			fprintf(f_synth," %4i",param.times);
			fprintf(f_synth," %.0f",zzz); // Mean number of evaluations

			zzz=0;

			for (i=0;i<param.times;i++)

			{

				zzz=zzz+(eval_run[i]-mean_eval)*(eval_run[i]-mean_eval);

			}

			zzz=sqrt (zzz/param.times);



			fprintf(f_synth," %f",zzz); // Standard deviation




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

			fprintf(f_synth, " %.0f",min_eps);
			
			printf("\nMEAN NUMBER OF EVALUATIONS: %.0f",mean_eval); // Mean number of evaluations
			
			printf(" +- %f",zzz); // Standard deviation
			printf(" [%.0f %.0f]",nb_eval_min,nb_eval_max);

			printf(" Min eps: %f",min_eps);
			printf("\nTOTAL NUMBER OF FAILURES: %i", failure_tot[level]);
			if (eval_f_sup>0) printf(" (after %.0f evaluations)",eval_f_sup);

			duration = duration_tot/(float)param.times;
			printf(" / MEAN TIME %4.4f",duration);

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

			if(model.guided_rand>0 && param.printlevel>0)
				display_density(param.DD,density_step);	
			
			if (param.printlevel>0)
			{
				printf("\n Evolution summary (eval. nb, best error value)");
				for (i=0;i<Max_histo;i++)
					printf("\n %.0f %f",best_eval[Max_histo-i-1][0][level],best_eval[Max_histo-i-1][1][level]);
			}
	
		}
		
		batch_runs=batch_runs-1;
		if (batch_runs>0) goto read_param;
		
		if (param.funct==12)  // Special case Coloring problem
		{
			fprintf(f_init_w,"\n-1"); // End of file
			fclose(f_init_w); // An init file  has been written.
		} 
		
		fprintf(f_visu0,"%i\n",-1); // End of visu file level 0
		fprintf(f_visu0,"\n %.0f evaluations \n",zzz);
		
		fprintf(f_visu1,"%i\n",-1);
		fprintf(f_swarm,"%i\n",-1); // End of save swarm file
		
return zzz;

}
// ----------------------------------------------------------------------------- ADAPT_ALPHA
float adapt_alpha(float m,struct model model)
{
float 	alpha,alpha_inf,alpha_sup;
float	k;
float lambda;
float	mu;

//return alea_float(model.alpha_min,model.alpha_max); // Random, for test purpose

alpha_inf=model.alpha_min;
alpha_sup=model.alpha_max;

if (alpha_sup<=alpha_inf) return alpha_inf;

alpha=(alpha_inf+alpha_sup)/2 + m*(alpha_sup-alpha_inf)/2;
return alpha;

/*
k=1;
mu=(alpha_inf+alpha_sup)/2;
lambda=(float)(mu-(alpha_sup-alpha_inf)/(exp(k)+1));
lambda=(alpha_sup-mu)/lambda;
alpha=(float)(mu+lambda*(alpha_inf+(1-1/(exp(k*m)+1))*(alpha_sup-alpha_inf) - mu));
return alpha;
*/
}

// ----------------------------------------------------------------------------- ADAPT_B
float adapt_b(float m, float alpha, struct model model)
{
float	b;
float	b_inf;
float	b_sup;
float	k;
float	lambda;
float	mu;

if (model.b_max<=model.b_min) return model.b_min; // Constant value


b_sup=MIN(model.b_max,(float)(alpha+1+2*sqrt(alpha))); // Max value for sure convergence
b_inf=MAX(model.b_min,(float)(alpha+1-2*sqrt(alpha)));

if (model.freedom==2) return alea_float(b_inf,b_sup);

b=(b_inf+b_sup)/2+m*(b_sup-b_inf)/2; 
//printf("\n m %f,  b %f",m,b);

//return b;

return alea_float(b_inf,b); // Partly random, for test purpose

// Partly random, for test purpose
b_inf=MAX(b_inf,0.8*b);
b_sup=MIN(1.2*b,b_sup);
//b=alea_float(b_inf,b_sup); 


return b;

/*
k=1;
mu=(b_inf+b_sup)/2;
lambda=(float)(mu-(b_sup-b_inf)/(exp(k)+1));
lambda=(b_sup-mu)/lambda;
b=(float)(mu+lambda*(b_inf+(1-1/(exp(k*m)+1))*(b_sup-b_inf) - mu));

//printf(" b %f",b);
return b;
*/
}

// ----------------------------------------------------------------------------- ADAPT_HOOD
float adapt_hood(int sw_size, float hood_size, float dhood, int type,int hood_min,int hood_max)
{
//float 	delta;
float 	new_hood_size;
float	mu=0.5;
float	N;

if (sw_size==hood_min) return (float)hood_min;
if (hood_max<=hood_min) return (float)hood_min;

N=(float)sw_size;

if (type==1)// Increase hood size
{
	/*
	delta=dhood*(1-(int)hood_size/(float)sw_size);
	new_hood_size=hood_size+delta;
	new_hood_size=MIN(new_hood_size,(float)hood_max);
	*/
	new_hood_size=hood_size+(N-hood_size)/(N-1);
}		
else // Decrease hood size
{
	/*
	delta=dhood*(int)hood_size/(float)sw_size;
	new_hood_size=hood_size-delta;
	new_hood_size=MAX(new_hood_size,(float)hood_min);
	*/
	new_hood_size=hood_size-(hood_size-1)/(N-1);;
}

if ((int)new_hood_size>Max_hood_size) Max_hood_size=(int)new_hood_size; // Max hood size, Just for info

//printf("\n type %i,sw size %i, hood %f, new hood %f",type, sw_size, hood_size,new_hood_size);

return new_hood_size;
}					
// ----------------------------------------------------------------------------- ADD_PARTICLE
struct swarm add_particle(struct swarm sw,struct param param, struct model model,int part_rank)
{
/*
Generate a new particle
part_rank is the rank of the particle which launched this generation
The bigger the swarm the less you can easily add a particle

*/
int				i;
struct swarm 	swt;

//if (E_p+E_k>=E0) return sw;
if (sw.size>=model.max_size) return sw;



n_add=n_add+1;
swt=sw;

i=sw.size;
//printf("\n Add a particle at rank %i, from particle %i",	i, part_rank);

if (model.gener_type==0 || param.funct!=12) // Normal case (not coloring problem)
{
	swt.part[i]=init_particle(param,1, model,sw,part_rank); // New particle
}

if (model.gener_type==1 && param.funct==12) // Special case for coloring problem
{
	swt.part[i]=gener_particle_color(sw,part_rank,param, model); // New particle
}


E_p=E_p+swt.part[i].pos.ep;
E_k=E_k+swt.part[i].vel.ek;

swt.size=sw.size+1;

return swt;
}

// ----------------------------------------------------------------------------- ALEA 
int alea(int min, int max) /* Random integer number between min and max */
{
int		ir;
float 	r;

r=alea_float(0,1);

ir=(int)(min+r*(max+1-min));
if (ir>max) ir=max;
return ir;
}

// ----------------------------------------------------------------------------- ALEA_DENSITY
float	alea_density(int density_step, int d,float xmin,float xmax,int type)
{
/* 
Guided Random Generation
[e-mail] 2001-09-20

Generate a random number for dimension d, according to the density of the positions already evaluated
Note: 
density[density_step][number_of_dimensions] is a global variable.
     Keep trace of the number of evaluated positions, for each dimension and each interval
density_step is a global variable. Number of intervals for dimension d
d is the dimension currently considered
xmin, xmax : min and max possible values for each dimension d
type = 0  => more probable where density is low
type = 1  => more probable where density is high
*/

float	dpr;
float	dx;
int		k;
int	i,j;
float	pr[10];
float	pr_c[11];
float	r;
float	sigma1, sigma2;


if (density_step<=1) return alea_float(xmin,xmax);

if (density_step>=10) printf("\n ERROR alea_density. density_step too high");

// Look for the max value
sigma1=0;
for (i=0;i<density_step;i++)
{
	if (sigma1<=(float)density[i][d]) sigma1=(float)density[i][d];	
}

if (sigma1==0) return alea_float(xmin,xmax);

k=1; // You may increase it for stronger effect

// Attraction/Repulsion function. YOU MAY CHANGE IT
if (type==0)// Attraction function (decreasing function of the density)
{	
	sigma2=0;
	for (i=0;i<density_step;i++)
	{
		pr[i]=(float)pow (sigma1-density[i][d],k);
		sigma2=sigma2+pr[i];
	}
}
else // Repulsion function (increasing function of the density)
{
	sigma2=0;
	for (i=0;i<density_step;i++)
	{
		pr[i]=(float)pow(density[i][d],k);
		sigma2=sigma2+pr[i];
	}
}


if (sigma2<=0) return alea_float(xmin,xmax);


// Cumulated
pr_c[0]=0;
pr_c[1]=pr[0];

for (i=2;i<=density_step;i++) 
{
	pr_c[i]=pr[i-1]+pr_c[i-1];
}

// Find a random number
r=alea_float(0,sigma2); 

// Modify it according to the attraction/repulsion function
dx=(xmax-xmin)/(float)density_step;

for (i=0;i< density_step;i++)
{
	j=density_step-i-1;
	if (pr_c[j]<=r)
	{
		dpr=pr_c[j+1]-pr_c[j];
		if (dpr==0) goto next_i;
		dpr=(r-pr_c[j])/dpr;
		r=xmin + dx*(j + dpr);
		return r;
	}
next_i:;
}

printf("\n ERROR alea_density. r = %f, j = %i, pr_c[j] = %f, sigma2 = %f",r,j,pr_c[j],sigma2);

printf("\n infinite waiting loop");infinite:goto infinite;
}
/*------------------------------------------------------------------- ALEA_DIFF */
int alea_diff(int min,int max, int num) /* Generate randomly an integer number # num */
{
int	temp1,temp2;

if (num==min)
	{
	temp1=alea(min+1,max);
	return temp1;
	}
	
if (num==max)
	{
	temp1=alea(min,max-1);
	return temp1;
	}

temp1=alea(min,num-1);
temp2=alea(num+1,max);

if (alea(0,100)<50)
	{
	return temp1;
	}
else
	{
	return temp2;
	}
}


// ----------------------------------------------------------------------------- ALEA_FLOAT 
float alea_float(float min, float max) /* Random  number between min and max */
{
float 	r;

/*
Note: rand_type is a global variable (cf models.txt )
0  => Use the "random" numbers generator
1  => Use the file of "random" numbers
2  => Use chaos. Seems not as good as randomness, even pseudo
*/
switch (rand_type)
{
	case 1:
	// Using this, you have the same result on different machines
	// BUT, sometimes performances are bad, probably because 
	// you  reuse the same pseudorandom numbers (possibility of cycles)

	r=rand_list[n_rand];
	n_rand=n_rand+1;
	if (n_rand>=max_rand-1)
	{
		printf("\n (reset rand_list after %i random numbers read)",max_rand);
		n_rand=0;
	}
	break;
	
	case 2: // The simplest "deterministic" way. But seems quite bad
		retry:
		chaos=4*chaos*(1-chaos); // Global variable
		// Distribution is well known to be not uniform (parabolic)
		// we keep only intermediate values
		if (chaos<0.2) goto retry; 
		if (chaos>0.8) goto retry;	
		r=(chaos-0.2)/0.6;
		break;
		
	default: // The classical and apparently the best one
		r=RAND_MAX; // Normally, RAND_MAX = 32767 = 2^31-1 
		r=rand()/r;
		break;
	}
	
	r=min+r*(max-min);
	return r;
	
}

// ----------------------------------------------------------------------------- ALL_DIFFERENT
struct position	all_different(struct position pos,struct param param)
{ 
/*
 Modify the position so that all components are different
 Should be used only when granularity>O
 Useful for some combinatorial problems, like Knapsack or Fifty-fifty
*/

int				d,d1;
float			dx;
int				k;
float			max_loop;
float			nloop;
struct position	post;

post=pos;

for (d=1;d<pos.size;d++) // For each component ...
{
	if (param.H.dx[d]>0) {dx=param.H.dx[d];} else {dx=param.eps;}
	max_loop=1+(param.H.max[d]-param.H.min[d])/dx;
	k=1;
	nloop=0;
	loop1:
	//... if component # of all previous, OK, ...
	for (d1=0;d1<d;d1++)
	{
		if (post.x[d]==post.x[d1]) // ... else, modify slightly and retry
		{
			loop2:
			post.x[d]=pos.x[d]-k*dx;
			if (k<0) {k=-k+1;} else {k=-k;}
			if (post.x[d]<param.H.min[d]) goto loop2;
			if (post.x[d]>param.H.max[d]) goto loop2;
				nloop=nloop+1;
				if (nloop==max_loop)