tsp.c

用C编写的粒子群寻优的优化算法

/***************************************************************************
                          TSP.c  -  description
                             -------------------
    begin                : mer mar 5 2003
    copyright            : (C) 2003 by 
    email                : [e-mail]
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
 //================================================= DISPLAY_GRAPH
void     display_graph(struct graph G)
{
  int i,j;
 for (i=0;i<G.N;i++)
 {
 printf( "\n%i /",i+1);
 for (j=0;j<G.N;j++)
    printf(" %.0f",G.v[i][j]);
 }
}

/*------------------------------------------------------------------- READ_GRAPH */
struct graph read_graph(FILE *file, int trace)
{
/* TSPLIB MATRIX format. One value/line. First value = number of nodes
each value = arc value. If <0 => no arc
*/
char			bidon[50];
char			comment[100];
double         delta;
char			edge_weight_format[30];
char			edge_weight_type[30];
struct graph	Gt;
int 			i,j;
char			name[20];
char			type[20];
float         zzz;


printf("\nI am reading the graph");

fscanf(file," %s %s\n",bidon,name);
fscanf(file," %s %s\n",bidon,type);
	fscanf(file," %s %s\n",bidon,comment);
	fscanf(file,"%s %i\n",bidon,&Gt.N); // dimension
	fscanf(file,"%s %s\n",bidon,edge_weight_type);
	fscanf(file,"%s %s\n",bidon,edge_weight_format);
	fscanf(file,"%s\n",bidon);


	printf("\n Name: %s, type: %s, (%s)",name,type,comment);
	printf("\n Number of nodes: %i",Gt.N);
	printf("\n %s %s\n",edge_weight_type,edge_weight_format);

	if (edge_weight_type[0]=='E' && edge_weight_format[0]=='F')
		{
		for (i=0;i<Gt.N;i++)
			{
                for (j=0;j<Gt.N;j++)
				{
                        fscanf(file,"%e ",&zzz);
                            Gt.v[i][j]=zzz;
					if (trace>2) printf(" %f",Gt.v[i][j]);
				}
			Gt.v[i][i]=0;
			}


    for (i=0;i<Gt.N;i++)
			{
			for (j=0;j<Gt.N;j++)
				{
                       if (Gt.v[i][j]>0) Gt.v[i][j]=integer_coeff*Gt.v[i][j];
                       }
                    }
  if (integer_coeff!=1) printf("\nWARNING. I have multiplied all arc values by %f",integer_coeff);

		return Gt;
		}
   printf("\nERROR by reading graph");
  return Gt;
}

//====================================================== SEQUENCE_SIMILAR
int   sequence_similar(struct seq s1,struct seq s2)
{
 /* Check if two sequences are globally similar:
 - same first value
 - same last value
 - same  set of values (not taking order into account)
 return 0 if false, 1 if true
 */
 int           i;
 int           j;
int            size;

   if (s1.s[0]!=s2.s[0]) return 0;

   size=s1.size;
   if (s1.s[size-1]!=s2.s[size-1]) return 0;
   if (size!=s2.size) return 0;

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

   for (i=1;i<size-1;i++)
   {
           for (j=0;j<size-1;j++)
           {
              if(s1.s[i]==s2.s[j]) goto next_i;
           }
            return 0;

    next_i:;
   }
  //printf("\n sequence_similar OK");
 return 1;
}

//====================================================== 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;

 }

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

//vt.size=alea(1,N);
vt.size=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,f;
int 			i,j,k;
int               improv;
int				levelling0;
int                  loop,loop_max;
struct particle	pt,pt1;
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;

 // loop_max=G.N;
 loop_max=1;
if (type==0) // Lazy Descent Method
	{
      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
			{
			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;
			}
		}
	return pt;
	}


 if (type==6) goto three_opt;

//same_best_thresh

if (type>=2 && type<=5) /*
				Check all immediate physical neighbours.
            If one is better, move,
            if (type=3 or 5) check again
			*/
	levelling0=levelling;
if (type>2) levelling0=1;

	v.size=1;
   improv=0;
	next_step2:
	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);

			if((levelling==0 && pt1.best_f<pt.best_f) || (levelling>0 && f_levelling(G,pt,pt1)<pt.fl)) // Move towards the best neighbour
				{
				if (trace>1)
					{printf("\n auto move => ");display_position(pt,1);}
				pt=pt1;
                     improv=1;
				if (type==3 || type==5) goto next_step2;
				}
			}
		}
      if (improv==0)
      {
      	if (trace>1) {printf("\n auto move => no move ");display_position(p,1);}
	return p;
      }
   else
	return pt;


 three_opt:
 if (trace>1)
    printf("\n auto_move, 3-opt");
 //scanf("%i",&b);
 v.size=3;
    for (b=0;b<G.N;b++)
    {
       v.comp[0][0]= p.x.s[b];

       for (d=0;d<G.N;d++)
       {
          v.comp[1][0]= p.x.s[d];
          v.comp[0][1]= p.x.s[d];

          for (f=0;f<G.N;f++)
          {
             v.comp[2][0]= p.x.s[f];
             v.comp[1][1]= p.x.s[f];

              v.comp[2][1]= p.x.s[b];
          pt=pos_plus_vel(pt,v,0,0,0);

		if(pt.f<p.f) // Stop as soon as a better position is found
		{
			return pt;
		}

           } // f
       } // d
    }  // b
    return pt;
}

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

	*/
{

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

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];
	for (j=i;j<GN;j++) // ... look for its rank in p1...
			{
			if (p1.x.s[j]==nt) // ... finds it at rank j

				{
				if (j!=i) // .. if it is not the same rank ...
					{
					pt.x.s[i]=n1; // ... change the node in pt ...
					vt.comp[vt.size][0]=n1;

								// ... then, looks 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;
							vt.comp[vt.size][1]=nt;
							vt.size=vt.size+1;

							if (vt.size>Vmax)
								goto stop;

							goto next_node;

							}
						}
					printf("\n ERROR. Don't find %i after rank %i in particle %i",n1,i+1, pt.rank);
					return vt;

					}
				}
			}
	next_node:;
	}

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;