particle_tools.c

Particle Swarm Optimization (PSO) for the TSP

//------------------------------------------------------------------- ALEA_VELOCITY 
struct velocity	alea_velocity(int N)
{
int				i;
struct velocity vt;

vt.size=alea(1,N);

for (i=0;i<vt.size;i++)
	{
	vt.comp[i][0]=alea(0,G.N-1);
	vt.comp[i][1]=alea_diff(0,G.N-1,vt.comp[i][0]);
	}
//display_velocity(vt);
return vt;
}

// ------------------------------------------------------------------- AUTO_MOVE 
struct particle auto_move(struct particle p,int type, int monotony, int levelling,int trace)
/* The particle moves slowly, looking for a better position */
{
int               b,d,g;
int               dd,gg;
double         best_before;
double         delta1,delta2,delta_max,delta_min;
double         eval_f0;
int 			i,j,k;
int               improv;
int				levelling0;
int                  loop,loop_max;
int               node;
struct particle	pt,pt0,pt1,pt2;
struct velocity v;

if (trace>1) 
	printf("\n Auto-move type %i for particle %i. Size %i \n",type,p.rank+1,p.x.size);
pt=p;
 if (pt.x.s[0]!=0) pt.x=rotate(pt.x,0);
 if (pt.best_x.s[0]!=0) pt.best_x=rotate(pt.best_x,0);

 // loop_max=G.N;
 loop_max=1;
if (type==0) // Lazy Descent Method  (random)
	{
      for(loop=0;loop<loop_max;loop++)
      {
	for (i=0;i<G.N;i++) // Move at random. 
		{
		pt.v=alea_velocity(1);
		pt=pos_plus_vel(pt,pt.v,monotony,1,levelling);

		if(pt.best_f<p.best_f) // Stop as soon as a better position is found
			{
			goto return_pt;
			}
		}
  }
	if (trace>1) {printf("\n auto move => no move ");display_position(p,1);}
	return p;	
	}
	
if (type==1) // Energetic
	{
	next_step1:
	for (i=0;i<G.N;i++) // Move at random
		{
		pt.v=alea_velocity(1);  // Just one transposition (switch)
		pt1=pos_plus_vel(pt,pt.v,monotony,1,levelling);

		if(pt1.best_f<pt.best_f) // Move as long as a better position is found
			{
			pt=pt1;
			goto next_step1;
			}
		}
	goto return_pt;	
	}	
	


 // type=2,  3, 5, 6, 7
 if (trace>1) printf("\n   transpositions");

/* 
				Check all immediate physical neighbours.
            If one is better, choose the best, and move
            if (type=3, 5, 6, 7 ) check again
*/

	v.size=1;
 
	next_step2:
    improv=0;
    delta_min=pt.best_f-target;
     
	for (j=0;j<G.N-1;j++)  // For each neighbour
	{
		v.comp[0][0]=j;
		for (k=j+1;k<G.N;k++) 
		{
			v.comp[0][1]=k;
			pt1=pos_plus_vel(pt,v,monotony,1,levelling);
				
               delta1=pt1.best_f-target;
               if( delta1<delta_min)
               {
                  if (trace>1)                     
					{printf("\n auto move => ");display_position(pt1,1);}
                  delta_min=delta1;
                  pt=pt1;  // Memorize the best new position
                  improv=1;
                  if (delta_min<=0) goto end_type2; // By any chance ...
			}
            }
	}

      if (improv==1 &&(type==3 || type==5 || type==6 || type==7)) // Move as long as a better position is found
      {
         goto next_step2;
      }   

             if (type==6) goto three_opt;
             if (type==7) goto three_opt_weak;
             
  end_type2:

	goto return_pt;

 //----------------------------------------------------------------------------------------------------------------
 three_opt:
   if(pt. best_time<splice_time ) goto return_pt;
 if (trace>1)
    printf("\n auto_move, 3-opt");
 
 three_opt_loop:
   // b=>g => d=> b, so  (b,d,g) => (g,d,b)
   // or
 // b=> d=>g=>b, so (b,d,g) => (d,g,b)
 
  delta_min=pt.best_f-target; // Maximum possible improvement

   // for (b=0;b<G.N;b++)
   for (b=0;b<G.N-2;b++)
    {
       //for (d=b+1;d<b+G.N;d++)
       for (d=b+1;d<G.N-1;d++)
       {
         //dd=d%G.N;
         //if (dd==b) continue;
          //for (g=dd+1;g<dd+G.N;g++)
          //for (g=dd+1;g<dd+G.N;g++)
          for (g=d+1;g<G.N;g++)
          {
         //gg=g%G.N;
         // if (gg==b) continue;
          //if (gg==dd) continue;

          pt1=pt; pt1.x=pt.best_x; // Try to directly modify the best previous position
          pt2=pt1;

          node=   pt1.x.s[b];
          pt1.x.s[b]=pt1.x.s[g];
          pt1.x.s[g]=pt1.x.s[d];
           pt1.x.s[d]=node;
           eval_f0=eval_f;
           pt1.f=f(pt1,-1,-1);
          eval_f=eval_f0+6.0/G.N; // This move could be done by just two transpositions
                                              // so it wouldn't be fair to count a complete evaluation                                   
           node= pt2.x.s[b];
           pt2.x.s[b]=pt2.x.s[d];
          pt2.x.s[d]=pt2.x.s[g];
           pt2.x.s[g]=node ;
           eval_f0=eval_f;
          pt2.f=f(pt2,-1,-1);
          eval_f=eval_f0+6.0/G.N;
          
                 delta1= pt1.f-target;
                  delta2= pt2.f-target;
                   if (MIN(delta1,delta2)>=delta_min) goto next_g;    // No improvement

                  if (delta1<delta2)
                  {
                  p1:
                  pt=pt1;
                  delta_min=delta1;
                  goto end_improve;
                  }

                  if (delta2<delta1)
                  {
                  p2:
                  pt=pt2;
                  delta_min=delta2;
                  goto end_improve;
                  }

                if (delta1==delta2)
                {
                  if (alea(0,1)==0) goto p1; else goto p2;
                }

           end_improve:

    next_g:
           } // g
       } // d
    }  // b

           if(pt.f<pt.best_f)
           {
              pt.best_x=pt.x;
             pt.best_f=pt.f;
              pt.best_time=time;
             if (move[4]>=3) BB_update(pt);
              goto  three_opt_loop;   // Try again
           }

    goto return_pt;


//-----------------------------------------------------------------------------------------------------
 three_opt_weak:    // Only three _consecutive_ nodes
  if(pt. best_time<splice_time ) goto return_pt;
  
 if (trace>1)
    printf("\n auto_move, 3-opt_weak");
    
 three_opt_weak_loop:

     delta_min=pt.best_f-target; // Maximum possible improvement
   // (b,b+1,b+2) => (b+2,b,b+1)
   // or
 // (b,b+1,b+2) => (b+1,b+2,b)
    for (b=0;b<G.N;b++)
    {
       dd=(b+1)%G.N;
       gg=(dd+1)%G.N;

       pt1=pt; pt1.x=pt.best_x; // Try to directly modify the best previous position
       pt2=pt1;
       
          node=   pt1.x.s[b];
          pt1.x.s[b]=pt1.x.s[gg];
          pt1.x.s[gg]=pt1.x.s[dd];
           pt1.x.s[dd]=node;
           eval_f0=eval_f;
           pt1.f=f(pt1,-1,-1);
          eval_f=eval_f0+6.0/G.N; // This move could be done by just two transpositions
                                              // so it wouldn't be fair to count a complete evaluation
  
          node=   pt2.x.s[b];
          pt2.x.s[b]=pt2.x.s[dd];
          pt2.x.s[dd]=pt2.x.s[gg];
           pt2.x.s[gg]=node ;
           eval_f0=eval_f;
          pt2.f=f(pt2,-1,-1);
          eval_f=eval_f0+6.0/G.N;

                  delta1= pt1.f-target;
                  delta2= pt2.f-target;
                   if (MIN(delta1,delta2)>=delta_min) goto next_b_2;    // No improvement
                   
                  if (delta1<delta2)
                  {
                  pt_1:
                  pt=pt1;
                  delta_min=delta1;
                  goto end_improve_2;
                  }
                  
                  if (delta2<delta1)
                  {
                  pt_2:
                  pt=pt2;
                  delta_min=delta2;
                  goto end_improve_2;
                  }                 

                if (delta1==delta2)
                {
                  if (alea(0,1)==0) goto pt_1; else goto pt_2;
                }
                
           end_improve_2:      
               if (trace>1)
                  printf("\n3-opt_weak =>%f",pt.f); 
  
    next_b_2:
		}  // next b
 
  if (pt.x.s[0]!=0) pt.x=rotate(pt.x,0);
     if(pt.f<pt.best_f)
     {
             pt.best_x=pt.x;
             pt.best_f=pt.f;
              pt.best_time=time;
             if (move[4]>=3) BB_update(pt);
 //             goto  three_opt_weak_loop;  // Try again
     }

  return_pt:
  //if (pt.x.s[0]!=0)
  pt.x=rotate(pt.x,0);
 //if (pt.best_x.s[0]!=0)
 pt.best_x=rotate(pt.best_x,0);
 pt.v=coeff_pos_minus_pos(pt,p,1,monotony,equiv_v,trace);
 return pt;

}

/*------------------------------------------------------------------- BEST_NEIGHBOUR */
struct particle	best_neighbour(struct swarm sw1,struct particle p,int k, int hood_type, int monotony,int equiv_v,int trace)
{
/* Find the k nearest particles of p, including p itself, in the swarm sw1, 
and find the best one according to the function value (f) 

*/

struct particle bestt;
double			dmin, dmax;
double			dist[Nmax];
int				i,j;
int				rank_dmin[Nmax];
double			kx;

if (hood_type==1) goto physical_hood;

// Case hood_type==0 (social hood)
/*--------------------- Neighbourhood. Super simple method: nearest by the rank */
/* Note: normally, this "social" hood should become a "physical" hood during the process */

kx=k;
if (best_hood_type==0)  // Best neighbour  in terms of current position
	{
	bestt=p;

	for (i=1;i<kx/2;i++)
		{
		j=p.rank+i;
		if (j>sw1.size-1) {j=j-sw1.size;};
	
		if(sw1.p[j].f<bestt.f)
			{
			bestt=sw1.p[j];
			}
		}
	for (i=1;i<kx/2;i++)
		{
		j=p.rank-i;
		if (j<0) {j=sw1.size+j;};
	
		if(sw1.p[j].f<bestt.f)
			{
			bestt=sw1.p[j];
			}
		}
	}


if (best_hood_type==1)   // Best neighbour in terms of best previous position
	{
	bestt=previous_best(p,0);

	for (i=1;i<kx/2;i++)
		{
		j=p.rank+i;
		if (j>sw1.size-1) {j=j-sw1.size;};
	
		if(sw1.p[j].best_f<bestt.f)
			{
			bestt=previous_best(sw1.p[j],0);
			}
		}
	for (i=1;i<kx/2;i++)
		{
		j=p.rank-i;
		if (j<0) {j=sw1.size+j;};
	
		if(sw1.p[j].f<bestt.f)
			{
			bestt=previous_best(sw1.p[j],0);
			}
		}
	}
	
if (trace>2)
	{
	printf("\n Best neighbour of %i is %i",p.rank+1,bestt.rank+1);

	}	
 //if (bestt.rank==p.rank) printf("\n Particle %i is its own best neighbour",p.rank);
 
return bestt;


/*--------------------- Physical hood */
physical_hood:

dmax=0;
	for (i=0;i<sw1.size;i++) // Compute distances p => particle i (including particle p itself
		{
		//if (i==p.rank) {dist[i]=1.2; continue;} // Note. distance can't be > 1
	
		if (i==p.rank) {dist[i]=0; continue;}
		dist[i]=distance(p,sw1.p[i],type_dist);
		if (dist[i]>dmax) dmax=dist[i];
		}	

	for (i=0;i<k;i++)
		{

		dmin=1.1*dmax;
		for (j=0;j<sw1.size;j++)  // Look for the k nearest (including p itself)
			{
			if (dist[j]<dmin)
				{
				dmin=dist[j];
				rank_dmin[i]=j;
				}                                                                                                         
			
			}
		dist[rank_dmin[i]]=1.1;
		}
		
if (best_hood_type==0)
	{	
	bestt=p;	
	for (i=0;i<k;i++) // Take the best of the k nearest (including p itself)
		{
		if (sw1.p[rank_dmin[i]].f<bestt.f)
			bestt=sw1.p[rank_dmin[i]];
		}
	}
	
if (best_hood_type==1)
	{	
	bestt=previous_best(p,0);	
	for (i=0;i<k;i++) // Take the best previous best of the k nearest (including p itself)
		{
		if (sw1.p[rank_dmin[i]].best_f<bestt.f)
			bestt=previous_best(sw1.p[rank_dmin[i]],0);
		}
	}	
	
return bestt;
}



/*------------------------------------------------------------- COEFF_POS_MINUS_POS */
struct velocity	coeff_pos_minus_pos(struct particle p1,struct particle p2,double coeff,int monotony,int equiv_v, int trace)
	/* Applying v to p2 gives p1, or, in short, p2+v=p1, or v=p1-p2 
	and, after coeff*v
	monotony = 1  => we have p1.f<=p2.f and we want, for each intermediate position p(t), 
							p(t).f<=p(t-1).f
	monotony = 2  => coeff*v is computed in order to keep a "monotony" according to the distance
						between the two particles

      equiv_v 	>0 = try to find a smaller equivalent velocity. 0 = do nothing. Cf equiv_vel()
	
	*/
{

int 			i,j,k,n1,nt;
int				GN;
struct particle pt;
struct velocity vt;
double			d0,d;

if (trace>2) printf("\n\n coeff_pos_minus_pos, particles %i and %i",p1.rank+1,p2.rank+1);

if (p1.rank>=0 && (p1.rank==p2.rank)) {vt.size=0; return vt;}; // If it is the same particle, of course the velocity is null

vt.size=0;
pt=p2;

GN=p2.x.size;

if (trace>2)
	{
	printf("\n coeff_pos_minus_pos, GN %i",GN);
	display_position(p1,0);
	display_position(p2,0);
	}


if (pt.x.s[0]!=0) pt.x=rotate(pt.x,0); /* To be sure the particle description begins on node 1 */
if (p1.x.s[0]!=0)  
	{
	printf("\n ERROR. coeff_pos_minus_pos. The first node of the particle %i is %i. Must be 1",p1.rank,1+p1.x.s[0]);
   return vt;
	}
   
 // Progressively update pt, in order to obtain p1
 
for (i=1;i<GN-1;i++) // For each node of pt, at rank i ...
	{
	nt=pt.x.s[i];
	n1=p1.x.s[i];
   if (n1==nt) continue; // It's OK
   
	for (j=i+1;j<GN;j++) // ... look for its rank in p1...
			{
			if (p1.x.s[j]==nt) // ... find it at rank j
				{
					pt.x.s[i]=n1; // ... change the node in pt ...
					vt.comp[vt.size][0]=n1; // Set the first part of the velocity component
					
								// ... then, look for the node of p1 in pt
					for (k=i+1;k<GN;k++)
					{
						if (pt.x.s[k]==n1)
							{
							pt.x.s[k]=nt;// Change the node
							vt.comp[vt.size][1]=nt;  // Set the second part of the velocity   component
							vt.size=vt.size+1;                                               

							if (vt.size>Vmax)
								goto stop;
								
							goto next_node;
					
                                    }
					} // next k
					printf("\n ERROR. coeff_pos_minus_pos. Don't find node %i after rank %i in particle %i",n1+1,i+1, pt.rank+1);
 
                          printf("\n particle pt"); display_position(pt,0) ;
                           printf("\n particle p2"); display_position(p2,0);
                            printf("\n particle p1"); display_position(p1,0);
                            display_velocity(vt);
                           return vt;
			} // end if (p1.x.s[j]==nt)
            }   // next j
	next_node:;
	} // next i
	
if (trace>2) {printf ("\n vt before coeff");display_velocity(vt);printf("\n");}
vt=coeff_times_vel(vt,coeff,equiv_v,trace);

//printf("\n monotony %i",monotony);

if (trace>2) {printf ("\n vt after coeff");display_velocity(vt);printf("\n");}

if (monotony!=2) 
	return vt;

//----------------------==========================================	
	d0=distance(p2,p1,type_dist); // Initial distance, with no coeff
if (d0==0) return vt;
if (coeff==1) return vt;

//printf("\n     dist. %f",d0);
if (coeff<1)
	{
	update_coeff:
	
	pt=pos_plus_vel(p2,vt,0,1,0);

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