particle_tools.c

Particle Swarm Optimization (PSO) for the TSP

shunt6: // "Spherical" distribution  around pg

   v2.size=distance(p,pg,type_dist);   // Radius
   v1.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(pi,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(pit,v1,0,1,0); // New position  pit+ c2*(pgt-pit)
   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, radius distance(pg,p) or distance(pi,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;

shunt10: // "Spherical" distribution  around pi, radius distance(pi,pg) or distance(pi,p)

  	v1.size=distance(pi,pg,type_dist);   // Radius
   v2.size=distance(pi,p,type_dist);
   v1.size=MAX(v1.size,v2.size);
      //v1.size=MIN(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(pi,v1,0,1,0);
     pt.v=coeff_pos_minus_pos(pt,p,1,0,equiv_v,trace);  // pt.v=pt-p
goto end;

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


 pt.x=rotate(pt.x,0); /* To be sure the tour description begins on node 1 */
 pt.best_x=rotate(pt.best_x,0);


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;
       pt.best_time=time;
	}
else
	{
	pt.best_x=p.best_x; // Don't change
	pt.best_f=p.best_f;
       pt.best_time=p.best_time;
	}

 pt.rank=p.rank;

if (tot_vel>0) effective_move=1; else effective_move=0;
BB_update(pt); // Update the blackboard
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
*/

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

int				last_trans;
int				n1,n2;
int               node;
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];
      found=0;
      for (i1=0;i1<G.N;i1++)
	{

		if (ptemp.x.s[i1]==n1)
         {
           rank1=i1;
           found=found+1;
           if (found==2) goto re_eval;
         }
 		if (ptemp.x.s[i1]==n2)
         {
           rank2=i1;
           found=found+1;
           if (found==2) goto re_eval;
         }                                                    
    }
  printf("\n ERROR. pos_plus_vel. Can't find node %i or node %i",n1, n2);        

  	
  re_eval:
  ptemp.x.s[rank1]=n2;   // was n1, becomes n2.
                                         // node no at rank1  _was_ at rank2
  ptemp.x.s[rank2]=n1;  // was n2, becomes n1
                                             // node no at rank2  _was_ at rank1
  
		if (type>0)// Just update the evaluation
			{
			ptemp.f=f(ptemp,rank1,rank2);  // nodes in ranks rank1 and rank2 have been switched
			
			if (step==1 && ptemp.f<ptemp.best_f)
				{
				ptemp.best_x=ptemp.x;// Memorize the best previous position
				ptemp.best_f=ptemp.f;
                      ptemp.best_time=time;
                      if (move[4]>=3) BB_update(ptemp);
				}
			if (monotony>0) ft[i+1]=ptemp.f;
			}
			
	}   // next i (velocity component)
		
	if (type==0)

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

			
		if (step==1 && ptemp.f<ptemp.best_f)
			{
			ptemp.best_x=ptemp.x;// Memorize the best previous position
			ptemp.best_f=ptemp.f;
                ptemp.best_time=time;
                if (move[4]>=3) BB_update(ptemp);
			}
		if (monotony>0) ft[i+1]=ptemp.f;
	
		}
//---------------------
//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:
pt.x=rotate(pt.x,0);
pt.best_x=rotate(pt.best_x,0);

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

}


/*------------------------------------------------------------------- 	PREVIOUS_BEST */
struct particle previous_best(struct particle p, int type)
{
struct particle pt;

if (type==0) // Generate a particle equal to the previous best
pt=p;
pt.x=p.best_x;

pt.f=p.best_f;
if (pt.x.s[0]!=0) pt.x=rotate(pt.x,0);
return pt;


// if (type==1)
// Memorize a previous best position in a particle structure

pt=p;
pt.best_x=p.x;
pt.best_f=p.f;
 pt.best_time=time;
return pt;
 
}

//====================================================== SEQUENCE_SUBST
struct   particle    sequence_subst(struct particle p1, struct particle p2,int size_max)
{/*
 s1= position (or best position) of p1, s2=best 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.
 subs_option is a globl 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;
 double  p1f;
 struct particle ptemp;
 struct  seq   s1,s2, ss1,ss2;
 int        sim;
 int        size;

 GN2=G.N/2;
 if (subst_option==0)
 {
    s1=p1.x;
    p1f=p1.f;
    }
  else
  {
     s1=p1. best_x;
     p1f=p1.best_f;
     if (p1.best_x.s[0]!=0) printf("\n WARNING. sequence_subst. Not node 0 at first place in best_x");
     if (p1.best_x.s[G.N]!=0) printf("\n WARNING. sequence_subst. Not node 0 at last place in best_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++) // Build a subsequence in s1, for particle p1
                                    // It begins on i in particle p1  
      {
           ij=(i+j)%G.N;
         ss1.s[j]=s1.s[ij];    
      }                                                                          

       for (k=0;k<G.N;k++)   // For the second particle ...
       {          
                  for (l=0;l<size;l++)   // ... build a subsequence in s2
                                                 // It begins on k in particle p2
                  {
                      kl=(k+l)%G.N;
                     ss2.s[l]=s2.s[kl];    
                  }
                  
                  sim=sequence_similar(ss1,ss2); // Compare the subsequences

                  if (sim==1)  // Globally the same
                  {    
                     // Check if it is better
                     ptemp.x=ss1;                                                        
                      f1=  f(ptemp,-1,-1);  // Warning. It  also adds the arc value ss1.s[size-1] => ss1.s[0]
                                                       // but it is not important, for it does exactly the same for f2.                  
                      ptemp.x=ss2;                 
                      f2=  f(ptemp,-1,-1);
   
                      if (f2<f1)  // Replace ss1 by ss2
                      {
                           for (j=0; j<size;j++)
                           {
                              ij=(i+j)%G.N;
                             if (subst_option==0) p1t.x.s[ij]=ss2.s[j];
                             else p1t.best_x.s[ij]=ss2.s[j];
                           }

                           if (subst_option==0) p1t.f=p1f-f1+f2;  // Update the f-value
                           else  p1t.best_f=p1f-f1+f2; 

                           goto end_subs; // It is not allowed to replace both the sequence
                                                   // _and_ its complementary: it would mean just replace
                                                   // p1.x (or p1.best_x) by p2.best_x. Not interesting.
                        }
      
              
                        if (p2.best_f-f2<p1f-f1) // The complementary susbsequence is better
                        {
                           if (subst_option==0) p1t.x=s2;
                           else   p1t.best_x=s2;
                           
                           for (l=0; l<size;l++)
                           {
                              kl=(k+l)%G.N;
                             if (subst_option==0) p1t.x.s[kl]=ss1.s[l]; // Keep ss1
                             else p1t.best_x.s[kl]=ss1.s[l]; 
                           }
                             if (subst_option==0) p1t.f=p2.best_f-f2+f1;   // Update the f-value
                             else  p1t.best_f=p2.best_f-f2+f1;
                           
                           goto end_subs;
                        }                 
                        goto end_sim;
                        
                    end_subs:
  
  
                    if (subst_option==0) p1t.x=rotate(p1t.x,0);
                     else p1t.best_x=rotate(p1t.best_x,0);
                     

                    // Stop as soon as there has been improvement
                      if (subst_option==0)
                      {
                         if (p1t.f<p1t.best_f)  // Memory
                         {
                               p1t.best_x=p1t.x;
                              p1t.best_f=p1t.f;
                               p1t.best_time=time;
                         }
                       }
                       
                         if(subst_option>0 && p1t.best_f<p1f)
                         {
                               p1t.best_time=time;
                         }
   
                         BB_update(p1t); // Update the blackboard
  
                            return p1t; // It has been improved
                       
                    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;
int   improv;
float  h2;
int   j;
int   k;
struct   particle pt;
int         size;

if (trace>1) printf("\n splice particle %i",p.rank);
  improv=0; 
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;  // Just a given rate of the swarm is checked
 
      for(k=0;k<size;k++)  // For each neighbour
      {
         if (k==p.rank) continue;
         //  If  comparing best_positions
         // and if none of the two best positions has changed since last splice_around call
         // => no need to check again.
         if(subst_option==1 && p. best_time<splice_time && sw.p[k].best_time<splice_time ) continue;
         
               pt=sequence_subst(p,sw.p[k],size_max);
               if (pt.f<p.f) improv=1;
               if(improv>0)
               {
                if (trace>1)
                   printf("  %f => %f (best= %f)",p.f,pt.f,pt.best_f);   
               return pt; // Stop asa there is some improvement
               }
          //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 of the previous best
      }
if (trace>1 && improv==0) printf("=> 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);
               if (pt.f<p.f) improv=1;
               if(improv>0)
               {
                if (trace>01)
                   printf(" \n %f => %f (best= %f)",p.f,pt.f,pt.best_f);
               return pt; // Stop asa there is some improvement
               }
      }
if (trace>1 && improv==0) printf("\n => 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;
}