tsp.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_coeff;

	}

if (coeff>1)
	{

	update_coeff2:

	pt=pos_plus_vel(p2,vt,0,1,0);
	d=distance(pt,p1,type_dist);


	if (d>=d0) // OK for monotony (increasing distance)
		{
		return vt;
		}
	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;
}

/*------------------------------------------------------------------- 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;
}


/*------------------------------------------------------------------- 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

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;

//------------------------------ 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*(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 (trace>2) printf("\n move_type %i",move_type);
/*
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;


shunt6: // "Spherical" distribution  around pg

  	v1.size=distance(p,pg,type_dist);   // Radius
   v2.size=distance(pi,pg,type_dist);
   v1.size=MAX(v1.size,v2.size);
if (v1.size==0)
{
   //printf("\n Null velocity for particle %i", p.rank) ;
   effective_move=0;
   return p;
}
   v1=alea_velocity(v1.size);
   pt=pos_plus_vel(pg,v1,0,1,0);
     pt.v=coeff_pos_minus_pos(pt,p,1,0,equiv_v,trace);  // pt.v=pt-p

goto end;

 shunt7: // "Spherical" distribution  around  (pi+pg)/2

     v1=coeff_pos_minus_pos(pi,pg,1,0,equiv_v,trace);  // v=pi-pg
       v2=coeff_times_vel(v1,0.5,equiv_v,trace);
   pt=pos_plus_vel(pg,v2,0,1,0);  //Center =pg+ (pi-pg)/2

   v1=alea_velocity(v1.size);
   pt=pos_plus_vel(pt,v1,0,1,0);
   pt.v=coeff_pos_minus_pos(pt,p,1,0,equiv_v,trace);  // pt.v=pt-p . New velocity (not really necessary)

goto end;

 shunt8: // "Spherical"  non uniform distribution  around  a point between pi and pg

v1.size=distance(pi,pg,type_dist);   // Radius

if (v1.size==0)
{
   //printf("\n Null velocity for particle %i", p.rank) ;
   effective_move=0;
   return p;

}

 c1=1/pi.f;
 c2=1/pg.f;
 c3=c1+c2;
 //c1=c1/c3;  // Weight for pi
 c2=c2/c3;  // Weight for pg

   v1=alea_velocity(v1.size);
   pit=pos_plus_vel(p,v1,0,1,0);      // A point around pi

   v1=alea_velocity(v1.size);
   pgt=pos_plus_vel(pg,v1,0,1,0);  // A point around pg

   v1=coeff_pos_minus_pos(pgt,pit,c2,0,equiv_v,trace); // c2*(pgt-pit)

   pt=pos_plus_vel(pi,v1,0,1,0); // New position
   pt.v=coeff_pos_minus_pos(pt,p,1,0,equiv_v,trace);  // New velocity (not really necessary)


goto end;


shunt9: // "Spherical" distribution  around p

  	v1.size=distance(p,pg,type_dist);   // Radius
   v2.size=distance(p,pi,type_dist);
   v1.size=MAX(v1.size,v2.size);
if (v1.size==0)
{
   //printf("\n Null velocity for particle %i", p.rank) ;
   effective_move=0;
   return p;
}
   v1=alea_velocity(v1.size);
   pt=pos_plus_vel(p,v1,0,1,0);
     pt.v=coeff_pos_minus_pos(pt,p,1,0,equiv_v,trace);  // pt.v=pt-p
goto end;


//--------------------------
end:


if (pt.x.s[0]!=0) pt.x=rotate(pt.x,0); /* To be sure the particle description begins on node 1 */


if (pt.f<p.best_f)  /* Save best previous position */
	{
	pt.best_x=pt.x; // Change to the current position
	pt.best_f=pt.f;
	}
else
	{
	pt.best_x=p.best_x; // Don't change
	pt.best_f=p.best_f;
	}

 pt.rank=p.rank;

if (tot_vel>0) effective_move=1; else effective_move=0;
return pt;
}

/*-------------------------------------------------------------------POS_PLUS_VEL */
struct particle	pos_plus_vel (struct particle p,struct velocity v, int monotony,int type,int levelling)
{
/* Move a particle according to its velocity
	particle + velocity => new particle

Modify: particle p
		eventually, extra_best
*/


double			ft[Vmax+1];
int				i,i1;

int				last_trans;
int				n1,n2;