particle_tools.c

Particle Swarm Optimization (PSO) for the TSP

	vt.size=vt.size-1;
//printf("\n  vt_size %i,   dist. %f",vt.size,d);	
	
	goto update_coeff2;
	}


//---------------
stop:	
/* Too many transpositions */

printf("\n *** WARNING. More than %i components (variable Vmax) for velocity of particle %i",Vmax, pt.rank);
printf("\n I have to 'cut' it");
return vt;


}

/*------------------------------------------------------------------- COEFF_TIMES_VEL */	
struct velocity	coeff_times_vel(struct velocity v,double coeff, int equiv_v,int trace)
{
double			coefft;
int				i;
int               new_vsize;
struct velocity vt,vt0;


if (v.size==0) return v;


vt=v;
if (coeff==0) {vt.size=0;return vt;};

coefft=coeff;

if (coeff<0)
	{
	// Inversion
	
	for (i=0;i<v.size;i++)
		{
		vt.comp[i][0]=v.comp[v.size-i-1][0];
		vt.comp[i][1]=v.comp[v.size-i-1][1];


		coefft=-coeff;
		}
	
	}

//printf("\n coeff %f",coefft);

 if (coefft==1) 
		return equiv_vel(v,equiv_v,trace);

if (coefft<1)
	{
	//vt.size=MAX(1,coefft*v.size);
	vt.size=coefft*v.size;
	return vt;
	}


//  coefft >1
vt0=vt;
for (i=0;i<coefft;i++)
	{
	vt=vel_plus_vel(vt,vt0,0,trace);
	}


 new_vsize=  coefft*v.size;
 if (new_vsize>Vmax)
 {
    vt.size=Vmax;
    printf("\nWARNING. Vmax is too small. Should be at least%i",new_vsize);
 }
	
if (equiv_v>0) vt=equiv_vel(vt,equiv_v,trace);

//printf("\n coeff_times_vel, vt.size %i",vt.size);

return vt;
}

/*------------------------------------------------------------------- 	COMPARE_PARTICLE */
int	compare_particle(struct particle p1,struct particle p2)
{
/*  p1 = p2  => 0
    p1 # p2  => 1  */


int	i;
int GN;
GN=p1.x.size;

for (i=0;i<GN;i++)
	{
	if (p1.x.s[i]!=p2.x.s[i]) return 1;
	}
	
if (p1.v.size!=p2.v.size) return 1;	
	
for (i=0;i<p1.v.size;i++)
	{
	if (p1.v.comp[i][0] != p2.v.comp[i][0]) return 1;

	if (p1.v.comp[i][1] != p2.v.comp[i][1]) return 1;
	}	
	
return 0;

}

/*------------------------------------------------------------------- DISTANCE */
double distance(struct particle p1,struct particle p2,int type_dist) // Distance function of common arcs
{
int	i,j;
int GN;
int n1,n2;
double sim;
struct velocity v;
double x;

GN=p1.x.size;


if (type_dist==0)  // "Distance" depending on the number of common arcs
	{
	sim=0;

	for (i=0;i<GN;i++)
		{
		n1=p1.x.s[i]; // Select an arc in p1
		n2=p1.x.s[i+1];
		for (j=0;j<GN;j++) // Check if it is in p2
			{
			if(p2.x.s[j]!=n1) goto next_node;
			if(p2.x.s[j+1]!=n2) goto next_node;
				sim=sim+1;
			next_node:;
			}
		}
	
	x=1-sim/GN;
	}
	
if (type_dist==1) // True distance 
	{
	v=coeff_pos_minus_pos(p1,p2,1,0,2,trace);
	x=v.size;
	}
	
return x;

}

/*------------------------------------------------------------------- EQUIV_VEL */
struct velocity equiv_vel(struct velocity v,int equiv_v,int trace)
{
/* Compute an equivalent and hopefully smaller velocity 
*/
struct particle p0,pt;
struct velocity vt;



if (equiv_v==0) return v;

// Reference point
if (equiv_v==1) p0=init;
if (equiv_v==2) p0=best_best;
//if (equiv_v==3) p0=extra_best;

p0.rank=-1;
pt=pos_plus_vel(p0,v,0,-1,0); // Move
vt=coeff_pos_minus_pos(pt,p0,1,0,0,trace); // Re-evaluate the velocity


if (trace>1)
	{
	if (vt.size!=v.size)
		{
		printf("\n equiv_vel. Initial v.size %i. Final v.size %i",v.size,vt.size);
		}
	}

if (vt.size<v.size)
	return vt;
else
	return v;
}




/*------------------------------------------------------------------- F_LEVELLING */
double	f_levelling(struct graph G, struct particle p1, struct particle p2)
{ /* Compute a temporary different function value, in order to 
  improve the convergence
p1 = particle to move
p2 = neighbour
*/

double			f_l;
int				i,j;
double			min1;
struct particle pt;
struct velocity vt;

if (p2.f<p1.f) return p2.f;

vt.size=1;
min1=p2.f;


for (i=0;i<G.N-1;i++)
	{
	for (j=i+1;j<G.N;j++)
		{
		vt.comp[0][0]=i; // Move 1 step
		vt.comp[0][1]=j;
		pt=pos_plus_vel(p2,vt,0,0,0); /* For each neighbour of p ...
								... looks for the best neighbour
								*/
		if (pt.f<min1) 
			min1=pt.f;

		}
	}


f_l=(min1+p1.f)/2;	 // Take the mean of the two bests

//printf("\n p1.f %f, p2.f %f, min1 %f, f_levelling %f",p1.f,p2.f,min1,f_l);

return f_l;
}



/*------------------------------------------------------------------- MOVE_TOWARDS */
struct particle	move_towards(struct particle p,struct particle pg,struct coeff c,int move_type,int explicit_vel,int conv_case,int equiv_v,int trace)
{
/*  

Basic equations for Particle Swarm Optimization (particular case with just one phi value)

		v(t+1) =alpha*v(t) +(beta*phi/2)*(pi-x) + (beta*phi/2)*(pg-x)			Equ. 1
		x(t+1) = x(t) + gamma*v(t) + ((1-(delta-eta*phi))/2)*(pi-x)
												((1-(delta-eta*phi))/2)*(pg-x)			Equ. 2
				
pi:  previous best position of the particle
pg:  global best position (best particle in the neighbourhoud)

or

		v(t+1) =alpha*v(t) +(beta*phi)*(pit-x)				Equ. 1
		x(t+1) = pit + gamma*v(t) + (eta*phi-delta)*(pit-x)	Equ. 2
pit=(pi+pg)/2

or

		v(t+1) =alpha*v(t) +(beta*phi)*(pit-x)						Equ. 1
		x(t+1) = x(t) + gamma*v(t) + (1-(delta-eta*phi))*(pit-x)	Equ. 2		


*/

struct particle pi;
struct particle pit,pgt;
struct particle	pt;

double         c1,c2,c3;
double			d1,d2;
double			alpha, beta, gamma, delta, eta;
double			phi;
int				tot_vel;
struct velocity	v1,v2,v3;
double			z0,z1;


if (trace>1) 
	printf("\n\n move particle %i",p.rank+1);

if (trace>2)
	display_position(p,1);

tot_vel=0;

pi=previous_best(p,0); // Generate the previous best particle
if (p.x.s[0]!=0) p.x=rotate( p.x,0);

alpha=c.c[0];
//alpha=alea_float(0,c.c[0]);
beta=c.c[1];
gamma=c.c[2];
delta=c.c[3];
eta=c.c[4];
//phi=alea_float(0,c.c[5]);
phi=c.c[5];
z0=beta*phi;
z1=(1-(delta-eta*phi));

if (conv_case==5) delta=(1-alpha)*phi;
if (conv_case==6) delta=1+(1-alpha)*phi;

pt.rank=p.rank;
//printf("\n coeffs %f %f %f %f %f %f",alpha, beta,gamma,delta,eta,phi);

if (explicit_vel>0) goto explicit_velocities;

v1=coeff_pos_minus_pos(pg,p,beta*phi,monotony,equiv_v,0);
pt=pos_plus_vel(p,v1,0,1,0);
tot_vel=tot_vel+v1.size;

goto end;
//----------------------
explicit_velocities:

if (move_type==1) goto shunt1;
if (move_type==2) goto shunt2;
if (move_type==3) goto shunt3;
if (move_type==4) goto shunt4;
if (move_type==5) goto shunt5;
if (move_type==6) goto shunt6;
if (move_type==7) goto shunt7;
if (move_type==8) goto shunt8;
if (move_type==9) goto shunt9;
if (move_type==10) goto shunt10;

//------------------------------ move_type=0

pt.v=coeff_times_vel(p.v,alpha,equiv_v,trace);
if (trace>2) {display_velocity(pt.v);printf(" alpha*v");}


v1=coeff_pos_minus_pos(pi,p,z0/2,monotony,equiv_v,trace);//	(beta*phi/2)*(pi-x)
if (trace>2) {display_velocity(pt.v);printf("(beta*phi/2)*(pi-x)");}

pt.v=vel_plus_vel(pt.v,v1,equiv_v,trace);
 
v1=coeff_pos_minus_pos(pg,p,z0/2,monotony,equiv_v,trace);//	(beta*phi/2)*(pg-x)
if (trace>2) {display_velocity(v1);printf(" (beta*phi/2)*(pg-x)");}

pt.v=vel_plus_vel(pt.v,v1,equiv_v,trace);//  New VELOCITY


//-----------


v1=coeff_times_vel(p.v,gamma,equiv_v,trace);// gamma*v
pt=pos_plus_vel(p,v1,0,1,0);


v1=coeff_pos_minus_pos(pi,p,z1/2,monotony,equiv_v,trace);//z1/2*(pi-x)
pt=pos_plus_vel(pt,v1,0,1,0);


v1=coeff_pos_minus_pos(pg,p,z1/2,monotony,equiv_v,trace);//z1/2*(pg-x)
pt=pos_plus_vel(pt,v1,0,1,0); // new POSITION

goto end;

//-------------------------
shunt1:

pt.v=coeff_times_vel(p.v,alpha,equiv_v,trace);
if (trace>2) {display_velocity(pt.v);printf(" alpha*v");}
pt=pos_plus_vel(p,pt.v,0,1,0);
tot_vel=tot_vel+pt.v.size;

pit=pos_plus_vel(pi,pt.v,0,1,0);
pgt=pos_plus_vel(pg,pt.v,0,1,0);


v1=coeff_pos_minus_pos(pit,pt,z0/2,monotony,equiv_v,trace);//	(beta*phi/2)*(pi-x)
if (trace>2) {display_velocity(pt.v);printf("(beta*phi/2)*(pi-x)");}

pt.v=vel_plus_vel(pt.v,v1,equiv_v,trace);
pt=pos_plus_vel(pt,v1,0,1,0);
tot_vel=tot_vel+v1.size;

pgt=pos_plus_vel(pgt,v1,0,1,0);

pgt.rank=Max_size;
v1=coeff_pos_minus_pos(pgt,pt,z0/2,monotony,equiv_v,trace);  //	(beta*phi/2)*(pg-x)
if (trace>2) {display_velocity(v1);printf(" (beta*phi/2)*(pg-x)");}

pt.v=vel_plus_vel(pt.v,v1,equiv_v,trace);//  New VELOCITY


//-----------


v1=coeff_times_vel(p.v,gamma,equiv_v,trace);// gamma*v
pt=pos_plus_vel(p,v1,0,1,0);
tot_vel=tot_vel+v1.size;

pit=pos_plus_vel(pi,v1,0,1,0);
pgt=pos_plus_vel(pg,v1,0,1,0);

pit.rank=Max_size;
v1=coeff_pos_minus_pos(pit,pt,z1/2,monotony,equiv_v,trace); //z1/2*(pi-x)
pt=pos_plus_vel(pt,v1,0,1,0);
tot_vel=tot_vel+v1.size;

pgt=pos_plus_vel(pgt,v1,0,1,0);
pgt.rank=Max_size;
v1=coeff_pos_minus_pos(pgt,pt,z1/2,monotony,equiv_v,trace); //z1/2*(pg-x)
pt=pos_plus_vel(pt,v1,0,1,0); // new POSITION
tot_vel=tot_vel+v1.size;


goto end;

//-----------------------
shunt2:
pt.v=coeff_times_vel(p.v,alpha,equiv_v,trace);
if (trace>2) {display_velocity(pt.v);printf(" alpha*v");}

v1=coeff_pos_minus_pos(pg,pi,0.5,monotony,equiv_v,trace); // (pg-pi)/2
pit=pos_plus_vel(pi,v1,0,1,0); // (pg+pi)/2;
tot_vel=tot_vel+v1.size;

pit.rank=Max_size;
v1=coeff_pos_minus_pos(pit,p,beta*phi,monotony,equiv_v,trace);//	(beta*phi)*((pg+pi)/2 - x)
if (trace>2) {display_velocity(pt.v);printf("(beta*phi)*((pg+pi)/2 - x)");}
 
v2=vel_plus_vel(pt.v,v1,equiv_v,trace);//  New VELOCITY

//-----------


v1=coeff_times_vel(p.v,gamma,equiv_v,trace);// gamma*v
pt=pos_plus_vel(pit,v1,0,1,0); // (pg+pi)/2 + gamma*v
tot_vel=tot_vel+v1.size;


v1=coeff_pos_minus_pos(pit,p,eta*phi-delta,monotony,equiv_v,trace);//(eta*phi-delta)*((pg+pi)/2 - x)

pt=pos_plus_vel(pt,v1,0,1,0); // new POSITION
pt.v=v2;
tot_vel=tot_vel+v1.size;

goto end;


//----------------------------------------------------------------
shunt3:  // This one is the nearest one to  "theoretical" PSO
/*
if (pg.rank<0) printf("\n pg.rank %i",pg.rank);
if (pi.rank<0) printf("\n pi.rank %i",pi.rank);
*/

v1=coeff_times_vel(p.v,alpha,equiv_v,trace);
if (trace>2) {display_velocity(v1);printf("%f*v",alpha);}

v3=coeff_pos_minus_pos(pg,pi,0.5,monotony,equiv_v,trace); // (pg-pi)/2

pit=pos_plus_vel(pi,v3,0,1,0); // (pg+pi)/2;
pit.rank=Max_size; // Arbitrary rank, so that coeff_pos_minus_pos doesn't give 0

v3=coeff_pos_minus_pos(pit,p,beta*phi,monotony,equiv_v,trace);//	(beta*phi)*((pg+pi)/2 - x)
if (trace>2) {display_velocity(v3);printf("(beta*phi)*((pg+pi)/2 - x)");}
 
v2=vel_plus_vel(v1,v3,equiv_v,trace);//  New velocity

if (trace>2) {display_velocity(v2);printf("v(t+1)");}
//-----------

if (conv_case==5) /* (Cf convergence_case). In this case, x(t+1)=(pg+pi)/2 + v(t+1)
					Just to speed up a bit the process */
	{
	pt=pos_plus_vel(pit,v2,0,1,0); // new POSITION
	pt.rank=p.rank;
	pt.v=v2; // New VELOCITY
	tot_vel=tot_vel+v2.size;
	goto end;
	}
	
if (conv_case==6) // x(t+1) = x(t) + v(t+1) 
	{
	pt=pos_plus_vel(p,v2,0,1,0); // new POSITION
	pt.rank=p.rank;
	pt.v=v2; // New VELOCITY
	tot_vel=tot_vel+v2.size;
	
	goto end;
	//       Some tests
	
	d1=distance(p,pg,type_dist);
	d2=distance(pt,pg,type_dist);
	if (d2>d1)
		{
		printf("\n Warning. Distance to best hood increased. %.0f => %.0f",d1,d2);
		
		}
	
	if (v2.size==0)
	{
	printf("\n particle doesn't move. alpha= %f, beta= %f, phi= %f",alpha,beta,phi);
	printf("\n velocity size %i",p.v.size); display_position(p,0);
	printf("\nprevious best"); display_position(pi,0);
	printf("\nbest hood"); display_position(pg,0);
	printf("\n(pg+pi)/2");display_position(pit,0);
	//v2=coeff_pos_minus_pos(pit,p,1,monotony,equiv_v,3);
	printf("\n((pg+pi)/2 - x)");display_velocity(v2);
	printf("\n(beta*phi)*((pg+pi)/2 - x)");display_velocity(v3);
	printf("\n");
	}
	
	
	goto end;
	}

v1=coeff_times_vel(p.v,gamma,equiv_v,trace);// gamma*v

v3=coeff_pos_minus_pos(pit,p,eta*phi-delta,monotony,equiv_v,trace);//(eta*phi-delta)*((pg+pi)/2 - x)
v1=vel_plus_vel(v1,v3,equiv_v,trace); 

pt=pos_plus_vel(pit,v1,0,1,0); // new POSITION

pt.v=v2; // New VELOCITY
tot_vel=tot_vel+v2.size;


goto end;

//----------------------------
shunt4: /* This one is the nearest one to "classical" PSO
			*/
v1=coeff_times_vel(p.v,alpha,equiv_v,trace);
v3=coeff_pos_minus_pos(pg,pi,0.5,monotony,equiv_v,trace); // (pg-pi)/2
pit=pos_plus_vel(pi,v3,0,1,0); // (pg+pi)/2;
pit.rank=Max_size;
tot_vel=tot_vel+v3.size;
v3=coeff_pos_minus_pos(pit,p,beta*phi,monotony,equiv_v,trace);//	(beta*phi)*((pg+pi)/2 - x)

v2=vel_plus_vel(v1,v3,equiv_v,trace);//  New velocity


v1=coeff_times_vel(p.v,gamma,equiv_v,trace);// gamma*v
v3=coeff_pos_minus_pos(pit,p,z1,monotony,equiv_v,trace); //(1-(delta-eta*phi))*(pit-x)
v3=vel_plus_vel(v1,v3,equiv_v,trace);
pt=pos_plus_vel(p,v3,0,1,0); // x(t+1)  new POSITION
tot_vel=tot_vel+v3.size;

pt.v=v2; // New VELOCITY

goto end;

// ---------------------------------
shunt5:  // Move towards the best of pi or pg

v1=coeff_times_vel(p.v,alpha,equiv_v,trace);

if (pi.f<pg.f) pit=pi; else pit=pg;


v3=coeff_pos_minus_pos(pit,p,beta*phi,monotony,equiv_v,trace);

v2=vel_plus_vel(v1,v3,equiv_v,trace);//  New velocity

//-----------

if (conv_case==5) /* (Cf convergence_case). I this case, x(t+1)=best_of_(pi,pg) + v(t+1)
					Just to speed up a bit the process */
	{
	pt=pos_plus_vel(p,v2,0,1,0); // new POSITION
	pt.v=v2; // New VELOCITY
	tot_vel=tot_vel+v2.size;
	goto end;
	}

v1=coeff_times_vel(p.v,gamma,equiv_v,trace);// gamma*v

v3=coeff_pos_minus_pos(pit,p,eta*phi-delta,monotony,equiv_v,trace);

pt=pos_plus_vel(pit,v1,0,1,0); // new POSITION
tot_vel=tot_vel+v1.size;

pt.v=v2; // New VELOCITY

goto end;