tsp.c

Particle Swarm Optimization (PSO) for the TSP

struct particle	pt,ptemp;
int				rank1,rank2;

struct velocity	vtemp;
// printf("\npos_plus_vel, v.size %i",v.size);

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

ptemp=p;
if (type>=0)
	{
	ft[0]=p.f;
	}

if (trace>2)
	{
	printf("\n pos_plus_vel, particle before:");
	if (type>=0) display_position(ptemp,1); else display_position(ptemp,0);
	display_velocity(v);
	}


for (i=0;i<v.size;i++)
	{
	n1=v.comp[i][0]; // Nodes to swap
	n2=v.comp[i][1];
	for (i1=0;i1<G.N+1;i1++)
		{

		if (ptemp.x.s[i1]==n1)
			{
			rank1=i1;
			if (rank1==G.N) rank1=0;
			ptemp.x.s[i1]=n2;
			}

		else
			{
			if (ptemp.x.s[i1]==n2)
				{
				rank2=i1;
				if (rank2==G.N) rank2=0;
				ptemp.x.s[i1]=n1;
				}
			}
		}


		if (type>0)// Just modif
			{
			ptemp.f=f(ptemp,rank1,rank2);

			if (step==1 && ptemp.f<ptemp.best_f)
				{
				ptemp.best_x=ptemp.x;// Memorize the best previous position
				ptemp.best_f=ptemp.f;
				}
			if (monotony>0) ft[i+1]=ptemp.f;
 /*
			if(ptemp.f<extra_best.f)
				{
				extra_best=ptemp;
				}
   */
			}

	}

	if (type==0)

		{
		ptemp.f=f(ptemp,-1,-1); // Complete


		if (step==1 && ptemp.f<ptemp.best_f)
			{
			ptemp.best_x=ptemp.x;// Memorize the best previous position
			ptemp.best_f=ptemp.f;
			}
		if (monotony>0) ft[i+1]=ptemp.f;
  /*
		if(ptemp.f<extra_best.f)
			{
			extra_best=ptemp;
			}
   */
		}
//---------------------
//printf("\npos_plus_vel,%f %f",p.best_f,ptemp.best_f);

if (monotony==0)
	{
	pt=ptemp;
	goto end;
	}
//----------------------
if (monotony==1) /* Slow down to keep the monotony on f value
					Note that monotony=1 supposes f(p)>=f(p+v)
					*/
	{
	pt=p;
	last_trans=v.size;
	vtemp.size=0;
	back:
	if (ft[last_trans]>ft[0])
		{

		if (trace>3) printf("\n %f > %f => back",ft[last_trans],ft[0]);
		last_trans=last_trans-1;
		vtemp.comp[vtemp.size][0]=v.comp[last_trans][0];
		vtemp.comp[vtemp.size][1]=v.comp[last_trans][1];
		vtemp.size=vtemp.size+1;
		goto back;
		}
		if (trace>3) display_velocity(vtemp);
		pt=pos_plus_vel(ptemp,vtemp,0,1,levelling); // Back
		goto end;
	}




printf("\n  ERROR in pos_plus_vel. Wrong value %i for monotony",monotony);
return p;

//---------------
end:
if (pt.x.s[0]!=0) pt.x=rotate(pt.x,0);
//if (extra_best.x.s[0]!=0) extra_best.x=rotate(extra_best.x,0);

if (trace>2)
	{
	printf("\n            particle after:");
	display_position(pt,1);
	}
return pt;
}


//====================================================== SEQUENCE_SUBST
struct   particle    sequence_subst(struct particle p1, struct particle p2,int size_max)
{/*
 s1= position of p1, s2=position of p2

 In s2, look for a sub sequence ss2 (size>=4) that is globally similar to one (ss1) in s1
 and so that s1/ss1=>ss2 is a better sequence (smaller total  length regarding graph G).
  f1 is the total length of the sequence s1.

 return: the new better sequence, after substitution

 G is a global variable.
 WARNING: each sequence is supposed to be a permutation of the G nodes, beginning in 0
 plus again the node 0 at the end (cycle).
 */
  double f1,f2;
  float  GN2;
 int        i;
 int        ij;
 int        im;
 int        j;
 int        k;
 int        km;
 int        kl;
 int        l;
 int        m;
 struct  particle p1t;
 struct particle ptemp;
 struct  seq   s1,s2, ss1,ss2;
 int        sim;
 int        size;
 double  zzz;

 GN2=G.N/2;
 s1=p1.x;
 s2=p2.best_x;
   p1t=p1;
  size=4;

  loop:
    ss1.size=size;  ss2.size=size;

   for (i=0;i<G.N;i++)    // Try all subsequences of size "size" and G.N-size
   {
      for (j=0;j<size;j++)
      {
           ij=i+j;
         if (ij>=G.N) ij=ij-G.N;
         ss1.s[j]=s1.s[ij];    // Build a subsequence in s1, for particle p1
                                    // It begins on i in particle p1
      }

       for (k=0;k<G.N;k++)  // For the second particle ...
       {
                  for (l=0;l<size;l++)
                  {
                      kl=k+l;
                      if(kl>=G.N) kl= kl-G.N;
                     ss2.s[l]=s2.s[kl];    // ... build a subsequence in s2
                                                // It begins on k in particle p2
                  }

                  sim=sequence_similar(ss1,ss2); // Compare the subsequences

                  if (sim==1)  // Globally the same
                  {
  //printf("\n size %i  ",size);
                     // Check if it is better
                     ptemp.x=ss1;
                     ptemp.x.s[size]=ss1.s[size-1];   // WARNING. All diagonal terms of G must be equal to 0
                                                                 // Normally, it is done in read_graph
                      f1=  f(ptemp,-1,-1);

                      ptemp.x=ss2;
                      ptemp.x.s[size]=ss2.s[size-1];
                      f2=  f(ptemp,-1,-1);

                      if (f2<f1)  // Replace ss1 by ss2
                      {

                           p1t.x=s1;
                           for (m=0; m<ss2.size;m++)
                           {
                              im=i+m;
                              if(im>=G.N) im=im-G.N;
                              p1t.x.s[im]=ss2.s[m];
                           }
                            p1t.f=p1.f-f1+f2;  // Update the f-value
                              if (trace>1)
                                 printf("\nsequence_subst. size %i, %f  => %f. ==> %f",size,f1,f2, p1t.f);

                           goto end_subs; // It is not allowed to replace both the sequence
                                                   // and its complementary: it would mean just replace
                                                   // p1.x by p2.best_x. Not interesting.
                        }


                        if (p2.best_f-f2<p1.f-f1) // The complementary susbsequence is better
                        {

                           p1t.x=s2;
                           for (m=0; m<ss1.size;m++)
                           {
                              km=k+m;
                              if(km>=G.N) km=km-G.N;
                              p1t.x.s[km]=ss1.s[m]; // Keep ss1
                           }
                           p1t.f=p2.best_f-f2+f1;   // Update the f-value
                            if (trace>1)
                               printf("\nsequence_subst. size %i, %f  => %f, ==> %f",G.N-size,p1.f-f1,p2.best_f-f2,p1t.f);
                           goto end_subs;
                        }
                        goto end_sim;
                    end_subs:

                        // Eventually update memory and return
                           // if ( p1t.f< p1t.best_f)
                           if ( p1t.f<best_best.best_f)
                           {
                             if (trace>1)
                                 printf("\n        sequence_subst %f ==> %f",p1.f,p1t.f);
                              p1t.best_x=p1t.x;
                              p1t.best_f=p1t.f;
                              // Eventually rotate
                              if (p1t.x.s[0]!=0) p1t.x=rotate(p1t.x,0);

                              return p1t;
                           }

                    end_sim:;
                     }   // end sim==1
       } // next k

   }    // next i

  size=size+1;
   if (size<=GN2) goto loop;  // Try a longer subsequence
   else
   {
     if (trace>1)
         printf("\n       sequence_subst. size [4...%i] ==> no improvement", size-1);
      return p1; // No possible improvement
   }

}

  //===================================================================== SPLICE_AROUND
struct   particle   splice_around(struct swarm sw, struct particle p, int hood_type,int hood_size)
 {
   /* The particle looks at all its neighbours (or the whole swarm),
   and tries to improve itself, by substituying a sequence in its position.

   "hood" is a global variable (neighbourhood size))
   "size_max" is global variable
    */

struct seq hood_rank;
int   i;
float  h2;
int   j;
int   k;
struct   particle pt;
int         size;

if (splice_around_option==1) goto hood;
 // The whole swarm
 // Look, compare and eventually move
// if (splice_around_option==0 )
     size=sw.size;
 //else  size=(G.N*splice_around_option)/100;
      for(k=0;k<size;k++)
      {
         if (k==p.rank) continue;
               pt=sequence_subst(p,sw.p[k],size_max);
          if (pt.best_f<best_best.best_f) return pt; // Stop asa enough improvement
         //if (pt.best_f<p.best_f) return pt; // Stop asa there is improvement
      }
 printf("\n splice=> NO improvement");
 return pt;

  hood:
 // Case hood_type==0 (social hood)
/*--------------------- Neighbourhood. Super simple method: nearest by the rank */
     if (hood_type!=0)
    printf("\nWARNING. splice_around: physical neighbourhood not yet written. I use social one");
// Build the neighbourhood
 k=0;
 h2=(float)hood_size/2;
	for (i=1;i<h2;i++)
		{
		j=p.rank+i;
		if (j>sw.size-1) {j=j-sw.size;}
         hood_rank.s[k]=j;
         k=k+1;
		}
	for (i=1;i<h2;i++)
		{
		j=p.rank-i;
		if (j<0) {j=sw.size+j;}
         hood_rank.s[k]=j;
         k=k+1;
		}

      hood_rank.size=k;
 // printf("\n hood_size %i, hood_rank.size %i",hood_size,k);
           //for (k=0;k<hood_rank.size;k++) printf("\n  hood_rank.s[%i]  %i",k,hood_rank.s[k]);

 // Look, compare and eventually move
      for(k=0;k<hood_rank.size;k++)
      {
        pt=sequence_subst(p,sw.p[hood_rank.s[k]],size_max);
 //printf("\n splice_around. %f %f %f %f",p.f,p.best_f,pt.f,pt.best_f);
          if (pt.best_f<best_best.best_f) return pt;
          // if (pt.best_f<p.best_f) return pt; // Stop asa there is improvement
      }
    printf("\n splice=> NO improvement");
  return pt;

 }
/*------------------------------------------------------------------- 	VEL_PLUS_VEL */
struct velocity	vel_plus_vel(struct velocity v1,struct velocity v2, int equiv_v,int trace) /* v1 _then_ v2 */
{

int				i;
struct velocity vt;

vt=v1;

vt.size=MIN(v1.size+v2.size,Vmax);

for (i=v1.size;i<vt.size;i++)
	{
	vt.comp[i][0]=v2.comp[i-v1.size][0];
	vt.comp[i][1]=v2.comp[i-v1.size][1];
	}

if (equiv_v>0) vt=equiv_vel(vt,equiv_v,trace);
return vt;
}

/*------------------------------------------------------------------- F */
double f(struct particle p,int rank1,int rank2) // Function to minimize
// Objective function for Traveling Salesman Problem
{
double				d;
double				delta_f;
double				ft;
int 				i;
int					n,n0,n1,n2,ni;
int					no_arc;
double				val1,val2;

eval_f=eval_f+(double)p.x.size/G.N;  // Sometimes, it is just a partial evaluation of the objective function
                                                         // so it wouldn't be fair to systematically add 1
no_arc=0;


if (rank1>=0) goto modif;

ft=0;

if (trace>3) printf("\n\n Compute f for particle %i",p.rank+1);

for (i=0;i<p.x.size;i++) /* For each arc in the particle ... */
	{
	n=p.x.s[i];
	n0=p.x.s[i+1];

// TEST
if (n0>G.N || n>G.N)
	{
	printf("\nERROR arc %i=>%i",n+1,n0+1);
	printf(" value %f",G.v[n][n0]);
	}


	/* ----------------- Penalty for unexistent arcs */
	if (G.v[n][n0]<0) /* If it is not an arc in the graph ... */
		{
		delta_f=tax_no_arc(no_arc,G); /* ... you pay a tax */
		no_arc=no_arc+1;/* ... (depending of how many no_arc have already been count)... */

//printf("\n no arc. Tax %f",delta_f);
		}
	else  /* If it is a "valid" arc ... */
		{
		delta_f=G.v[n][n0]; /* ...you pay the right price ... */
		}
	ft=ft+delta_f;

	}
goto end;


//--------------------------
modif: /* IMPORTANT. This works only when taxes_noarc is a CONSTANT
			Also, it supposes G(i,i)=0
		 In the particle p, just two nodes, at ranks rank1 and rank2,
		   has been swapped. It is less costly to just modify the f value
		*/
n1=p.x.s[rank1];
n2=p.x.s[rank2];
ft=p.f;
//display_position(p,0,G);
//printf("\n f before %f, swapped nodes  %i, %i, ranks %i and %i",p.f,n1+1,n2+1, rank1,rank2);

if (rank1==0) i=G.N-1; else i=rank1-1;
ni=p.x.s[i];
// arc p(i)=>n1  was arc p(i)=>n2
if (G.v[ni][n1]<0) val1=tax_no_arc(no_arc,G); else val1=G.v[ni][n1];
if (G.v[ni][n2]<0) val2=tax_no_arc(no_arc,G); else val2=G.v[ni][n2];

ft=ft-val2+val1;


if (rank1==G.N-1) i=0; else i=rank1+1;
ni=p.x.s[i];
// arc n1=>p(i)  was arc n2=>p(i)
if (G.v[n1][ni]<0) val1=tax_no_arc(no_arc,G); else val1=G.v[n1][ni];
if (G.v[n2][ni]<0) val2=tax_no_arc(no_arc,G); else val2=G.v[n2][ni];

ft=ft-val2+val1;

//---------

if (rank2==0) i=G.N-1; else i=rank2-1;
ni=p.x.s[i];

	// arc p(i)=>n2  was arc p(i)=>n1
	if (G.v[ni][n1]<0) val1=tax_no_arc(no_arc,G); else val1=G.v[ni][n1];
	if (G.v[ni][n2]<0) val2=tax_no_arc(no_arc,G); else val2=G.v[ni][n2];

	ft=ft-val1+val2;


if (rank2==G.N-1) i=0; else i=rank2+1;
ni=p.x.s[i];

	// arc n2=>p(i)  was arc n1=>p(i)
	if (G.v[n1][ni]<0) val1=tax_no_arc(no_arc,G); else val1=G.v[n1][ni];
	if (G.v[n2][ni]<0) val2=tax_no_arc(no_arc,G); else val2=G.v[n2][ni];
	ft=ft-val1+val2;

if (abs(rank1-rank2)==1 || abs(rank1-rank2)==G.N-1) // Particular cases n1=>n2 or n2=>n1
	{
	if (G.v[n1][n2]<0) val1=tax_no_arc(no_arc,G); else val1=G.v[n1][n2];
	if (G.v[n2][n1]<0) val2=tax_no_arc(no_arc,G); else val2=G.v[n2][n1];
	ft=ft-val1-val2;
	}

//printf("   f after %f",ft);


end:

/*
if (ft<extra_best.f)
	{
	extra_best=p;
	extra_best.f=ft;
	if (extra_best.x.s[0]!=0) extra_best.x=rotate(extra_best.x,0);
	}
*/



return ft;

// Just for partial Map
d=distance(best_solution,p,type_dist);
	fprintf(f_trace,"%f %f\n",d,ft);

return ft;

}

/*------------------------------------------------------------------- 	TAX_NO_ARC */
double tax_no_arc(int l,struct graph G)
	{
	/*
	Penalize the l-th unexistent arc
	*/

double 	tax=1.1; /* >=1 */
double 	x;

	//x=tax*G.l_max+(G.N-l)*(G.l_max-G.l_min);
	x=tax*G.l_max+(G.N-1)*(G.l_max-G.l_min);

	return x;

	}

/*------------------------------------------------------------------- DISPLAY_POSITION */
void display_position(struct particle p,int type)
{
int 	i;
float 	x;

printf("\n particle %i: ",p.rank+1);

if (type>1) goto best_previous;

for (i=0;i<p.x.size+1;i++)
	printf(" %i",p.x.s[i]+1);

if (trace>1)
	{
	printf("\n arc values: ");
	for (i=0;i<p.x.size;i++)
		{
		x=G.v[p.x.s[i]][p.x.s[i+1]];
		if (x<0) printf("? "); else printf("%f ",x);
		}
	}

	 printf(" f value:%f",p.f);
	printf(" v.size: %i,",p.v.size);

   if (trace>2)
	display_velocity(p.v);

return;

best_previous:
	printf("\n best previous: ");
	for (i=0;i<p.best_x.size+1;i++)
		printf(" %i",p.best_x.s[i]+1);

   	 printf(" f value:%f",p.best_f);

}
/*------------------------------------------------------------------- DISPLAY_SWARM */
void display_swarm(struct swarm sw)
{
int i,i_best;

i_best=0;
printf("\n SWARM");
for (i=0;i<sw.size;i++) /* For each "particle"  */
	{
	display_position(sw.p[i],1);
	if (sw.p[i].f<sw.p[i_best].f) i_best=i;

	}
printf("\n Best particle: %i, f value %.1f",i_best+1,sw.p[i_best].f);
printf("\n");


}
/*------------------------------------------------------------------- DISPLAY_VELOCITY */
void display_velocity(struct velocity v)
{
int  i;

printf("\n  velocity (size=%i) :",v.size);

for (i=0;i<v.size;i++)
	{
	printf(" %i<=>%i ",v.comp[i][0]+1,v.comp[i][1]+1);
	}
}