pso8.c

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

				{
					printf("\n WARNING. Some components are equal (%i and %i) Impossible to modify this particule:",d1,d);
					display_position(post);
					return post;
				}
				goto loop1;
		}
	}
}
return post;
}

// ----------------------------------------------------------------------------- BEST
int best(struct swarm sw)// Find the "previous best" particle in the whole swarm
{
int i;
int	sw_best;

sw_best=0; // First, arbitrary value for the rank of the best "previous best"

for( i = 1;i<sw.size;i++)  // Look for a better one
	{
	if (sw.part[i].prev_best.f<sw.part[sw_best].prev_best.f)
		sw_best=i;
	}	
return sw_best;
}

// ----------------------------------------------------------------------------- BETTER_THAN
int	better_than(struct particle part1,struct particle part2,float eps)
{ // return 1 if part1 is better than part2 (or same quality)

if (part1.prev_best.f>part2.prev_best.f+eps) return 0; // pos2 is better
if (part1.prev_best.f<part2.prev_best.f-eps) return 2; // pos1 is better

if ((part1.	prev_prev_best_f-part1.prev_best.f)>=(part2.prev_prev_best_f-part2.prev_best.f)) return 1;

return 0;

}

// ----------------------------------------------------------------------------- COMPLETE_PART 
struct particle complete_part(struct particle par, struct param param, struct model model)	
 // WARNING: position and velocity must be already initialized. Complete the particle
{
int				d;
struct particle	part;
float			total;

part=par;

part.pos.label=label[level];
part.prev_best.label=label[level];
label[level]=label[level]+1;

part.pos.size=param.DD;
part.vel.size=param.DD;
part.hood_size=(float)model.hood_min;


if (param.printlevel>2)
printf("\n complete_part, label %.0f, H %i",part.pos.label, part.H);

part=constrain(part,param,model); // Keep in the search space

	
if (param.all_diff>0) // Eventually, modify so that all components are different
		part.pos=all_different(part.pos,param);
	
total=tot(part.pos, param, model);       //evaluation of the position
part.pos.f=eval(total,param.target);

part.pos.X1=X1; // For Leanordo visualization
part.pos.X2=X2;

part.prev_best=part.pos; // previous best = current position

for (d=0;d<Max_histo;d++) // previous f-values for the last  steps 	
{
	part.pos.prev_f[d]=part.pos.f;
}
	
part.prev_prev_best_f=part.pos.f;
part.pos.ep=energy_p(part,param.eps);	

if(part.pos.f<param.eps) 
	{
	best_result[level]=part.pos;
	chance[level]=1;
	}

part.vel.combin=param.combin; // Remember what kind of velocity is used (combinatorial or not)
part.vel.ek=energy_k(part,param.H); // Kinetic energy

return part;
}

// ----------------------------------------------------------------------------- CONSTRAIN
struct particle constrain(struct particle part,struct param param, struct model model)
{ /*
 Keep the particle in the definition domain 
*/
int				cmax;
int				d;
struct particle partt;

float			xmax;
float			xmin;	

//if (param.combin>0) return part;

cmax=(int)param.max_fr;
partt=part;

if (param.printlevel>2)
	printf("\n constrain.  label %.0f, H %i",partt.pos.label,partt.H);

if (param.combin==0) // If not combinatorial search
{
	partt.pos=domain_pos(partt.pos,param); // Take granularity into account	
}

if (param.printlevel>2)
	printf("\n Granularity taken into account");
	
if (param.funct==12) // Graph Coloring problem
{
	for (d=0;d<part.pos.size;d++) // Keep in the box
	{
		if( part.pos.x[d]>param.H.max[d])
		{ 
			partt.pos.x[d]=param.H.max[d];
			continue;
		}

		if( part.pos.x[d]<param.H.min[d])
		{ 
			partt.pos.x[d]=param.H.min[d];
		}
	}

	if (model.H==1) // Coloring/FAP. As few unsatisfied constraints as possible
	{
		xmin=param.H.min[0]; // Min admissible color
		xmax=param.H.max[0]; // Max admissible color
		if (model.confin<4) // Methods 0 to 3
		{
			partt=admissible_color(partt,M,model.confin,cmax); // (Note: matrix M is a global variable)
		}


		if (model.confin==4) // Method 4
		{
			partt=minimize_constraint(partt);
		}
		
		if (model.confin==5) // Method 5
		{
			partt.pos= color_6(M,dplus,part.pos,param.printlevel); // Protein folding
		}
	}
	
	if (model.H==2)  // As few colors as possible
	{
		if (model.confin>0)
		{
			partt=minimize_color(partt,cmax);
		}
	}

	if (model.H==3)
	{
		if (model.confin>0)
		{		
			if (partt.H==1) 
			{
				partt=minimize_constraint(partt);
				//partt=admissible_color(partt,M,model.confin,xmin,xmax,cmax);
			}

			if (partt.H==2) 
			{
				partt=minimize_color(partt,cmax);
			}
			
			partt.H=3-partt.H; // Just for Coloring problem. Switch the projection subspace
		}	
	}

return partt;
}


// Normal PSO	
for (d=0;d<partt.pos.size;d++)
{
	if( partt.pos.x[d]>param.H.max[d])
	{ 
		partt.pos.x[d]=param.H.max[d];
		partt.vel.v[d]=-part.vel.v[d]/2;
		//partt.vel.v[d]=0;
		continue;
	}
		
	if( partt.pos.x[d]<param.H.min[d])
	{ 
		partt.pos.x[d]=param.H.min[d];
		partt.vel.v[d]=-part.vel.v[d]/2;
		//partt.vel.v[d]=0;	
	}
}	
return partt;
}


// ----------------------------------------------------------------------------- DISPLAY_DENSITY
void display_density(int dim,int density_step)
{
// Note: density[][] is a global variable

int 	d;
int 	i;
float sigma;

printf("\n Density matrix. One line / dimension");

for (d=0;d<dim;d++)
{
	printf("\n");
	sigma=0;
	for (i=0;i<density_step;i++)
	{
		printf(" %3i",density[i][d]);
		sigma=sigma+density[i][d];
	}
	printf("  = %.0f",sigma);
}
}


// ----------------------------------------------------------------------------- DISPLAY_MODEL
void display_model(struct model model, int level)
{
printf("\n------\n Model level %i",level);
printf("\n alpha...............:[%1.3f,%1.3f] ",model.alpha_min,model.alpha_max);
printf("\n b ..................:[%1.3f,%1.3f] ",model.b_min,model.b_max);
printf("\n hood_size...........:[%i,%i] ",model.hood_min,model.hood_max);
printf("\n init_size...........: %i ",model.init_size);
printf("\n min_size............: %i ",model.min_size);
printf("\n selection...........: %i ",model.selection);
printf("\n gener_type..........: %i ",model.gener_type);
printf("\n init_type...........: %i ",model.init_type);
printf("\n H ..................: %i ",model.H);
printf("\n confinement ........: %i ",model.confin);
printf("\n parallel ...........: %i ",model.parallel);
printf("\n local_search .......: %i ",model.local_search);
printf("\n rehope .............: %i ",model.rehope);
printf("\n move_type...........: %i ",model.move_type);
printf("\n local queens........: %i ",model.queen);
printf("\n guided_rand.........: %i ",model.guided_rand);
printf("\n rand_type...........: %i ",model.rand_type);

if (model.freedom==1) printf("\n Before a move, all components are modified the same way");
if (model.freedom==2) printf("\n Before a move, each components is modified independently");

}

// ----------------------------------------------------------------------------- DISPLAY_PARAM
void display_param(struct param param)

{
int 	d;
float 	p1,p2,p3;
int 	times;

printf("\n\n PROBLEM");
printf("\n funct...........: %2i",param.funct);
printf("\n Dimension.......: %2i",param.DD);
printf("\n Search space:");
printf("\n xmin    xmax    dx (granul)");
p1=0;p2=0;
p3=-1;// Arbitrary unused value for dx

for (d=0;d<param.H.dim;d++)
	{
	if (param.H.min[d]!=p1 || param.H.max[d]!=p2 || param.H.dx[d]!=p3)
		{
		if (d>0) printf(" (%i times)",times);
		p1=param.H.min[d];
		p2=param.H.max[d];
		p3=param.H.dx[d];
		printf("\n %f  %f  %f",p1 ,p2,p3);
		times=1;
		}
	else
		{
		times=times+1;
		}
	}
printf(" (%i times)",times);
printf("\n target..........: %4.3f",param.target);
printf("\n eps.............: %4.10f",param.eps);
printf("\n printlevel......: %2i",param.printlevel);
printf("\n times...........: %2i",param.times);
printf("\n parapet.........: %f",param.parapet);
if (param.all_diff==1) printf("\nAll components of the solution must be different");
if (param.combin==1)
{
	 printf("\n I do use combinatorial kit");
}

else
{
	printf("\n I don't use combinatorial kit ");
}
printf("\n---------\n");
}

// ----------------------------------------------------------------------------- DISPLAY_POSITION
void display_position(struct position pos)
{
int				d;
int				ncar=0;
struct position post;

post=pos;

printf("\n position size %i\n ",pos.size);
for (d=0;d<pos.size;d++)
{
	printf("%6.6f ",post.x[d]);
	ncar=ncar+10;
	if (ncar>75)
	{
		ncar=0;
		printf("\n");
	}
}
printf (" abs(f-target) = %3.4f",post.f);

}

// ----------------------------------------------------------------------------- DISPLAY_SWARM
void display_swarm(struct swarm sw,int type)
{
int i;

if (type==1)
{
	printf("\n Current swarm at level %i", level);
	for (i=0;i<sw.size;i++)
	{
		display_position(sw.part[i].pos);
		printf("  best: %f",sw.part[i].prev_best.f);
	}
}
if (type==2)
{
	printf("\n Best previous swarm at level %i", level);	
	for (i=0;i<sw.size;i++)
	{
		display_position(sw.part[i].prev_best);
		printf("  new: %f",sw.part[i].pos.f);
	}
}
	
if (type==3)
{
	printf("\n Particle f values");
	for (i=0;i<sw.size;i++)
	{
		printf("\n %i  best: %f new: %f",i,sw.part[i].prev_best.f,sw.part[i].pos.f);
	}
}	

}

// ----------------------------------------------------------------------------- DISPLAY_VELOCITY
void display_velocity(struct velocity vel)
{

int d;

printf("\n velocity. combin = %i",vel.combin);
printf("\n  ");
for (d=0;d<vel.size;d++)
{
	printf("%3.2f ",vel.v[d]);
}


if (vel.combin==1) // "Combinatorial" velocity
{
	printf("\n");
	for (d=0;d<vel.size;d++)
	{
		if (vel.v2[d]==NA)
		{
			printf("    NA");
			continue;
		}
		if (vel.v2[d]==NO)
		{
			printf("    NO");
			continue;
		}
		if (vel.v2[d]==AS)
		{
			printf("    AS");
			continue;
		}
		
		printf("% 3.2f",vel.v[d]);
	}	
}


}

// ----------------------------------------------------------------------------- DISTANCE
float	distance(struct position pos1,struct position pos2)
{
float 	d,dx;
int		i;
d=0;
for (i=0;i<pos1.size;i++)
{
	dx=pos1.x[i]-pos2.x[i];
	d=d+dx*dx;
}
d=(float)sqrt(d);
return d;
}

// ----------------------------------------------------------------------------- DISTANCE_MIN
float	distance_min(struct swarm sw, struct position pos)
{
/*
Let's consider the  set of points S:
- all particle positions

We compute the minimum distance from pos to S

*/
float			dist;
float			dist_min;
int				i;

dist_min=(float)-NA; // Arbitrary too high value

for (i=0;i<sw.size;i++)
{
	dist=distance(pos,sw.part[i].pos);
	if (dist<dist_min) dist_min=dist;
}
return dist_min;
}

// ----------------------------------------------------------------------------- DIVERG