ls.c

蚁群算法解决TSP问题的源代码

/*

       AAAA    CCCC   OOOO   TTTTTT   SSSSS  PPPPP
      AA  AA  CC     OO  OO    TT    SS      PP  PP
      AAAAAA  CC     OO  OO    TT     SSSS   PPPPP
      AA  AA  CC     OO  OO    TT        SS  PP
      AA  AA   CCCC   OOOO     TT    SSSSS   PP

######################################################
##########    ACO algorithms for the TSP    ##########
######################################################

      Version: 1.0
      File:    ls.c
      Author:  Thomas Stuetzle
      Purpose: implementation of local search routines
      Check:   README and gpl.txt
      Copyright (C) 1999  Thomas Stuetzle
*/

/***************************************************************************

    Program's name: acotsp

    Ant Colony Optimization algorithms (AS, ACS, EAS, RAS, MMAS, BWAS) for the 
    symmetric TSP 

    Copyright (C) 2004  Thomas Stuetzle

    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.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    email: stuetzle no@spam informatik.tu-darmstadt.de
    mail address: Universitaet Darmstadt
                  Fachbereich Informatik
                  Hochschulstr. 10
                  D-64283 Darmstadt
		  Germany

***************************************************************************/

#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <limits.h>

#include "InOut.h"
#include "TSP.h"
#include "ants.h"
#include "utilities.h"

long int ls_flag;          /* indicates whether and which local search is used */ 
long int nn_ls;            /* maximal depth of nearest neighbour lists used in the 
			      local search */ 
long int dlb_flag = TRUE;  /* flag indicating whether don't look bits are used. I recommend 
			      to always use it if local search is applied */

long int * generate_random_permutation( long int n )
/*    
      FUNCTION:       generate a random permutation of the integers 0 .. n-1
      INPUT:          length of the array
      OUTPUT:         pointer to the random permutation
      (SIDE)EFFECTS:  the array holding the random permutation is allocated in this 
                      function. Don't forget to free again the memory!
      COMMENTS:       only needed by the local search procedures
*/
{
   long int  i, help, node, tot_assigned = 0;
   double    rnd;
   long int  *r;

   r = malloc(n * sizeof(int));  

   for ( i = 0 ; i < n; i++) 
     r[i] = i;

   for ( i = 0 ; i < n ; i++ ) {
     /* find (randomly) an index for a free unit */ 
     rnd  = ran01 ( &seed );
     node = (long int) (rnd  * (n - tot_assigned)); 
     assert( i + node < n );
     help = r[i];
     r[i] = r[i+node];
     r[i+node] = help;
     tot_assigned++;
   }
   return r;
}

void two_opt_first( long int *tour ) 
/*    
      FUNCTION:       2-opt a tour 
      INPUT:          pointer to the tour that undergoes local optimization
      OUTPUT:         none
      (SIDE)EFFECTS:  tour is 2-opt
      COMMENTS:       the neighbourhood is scanned in random order (this need 
                      not be the best possible choice). Concerning the speed-ups used 
		      here consult, for example, Chapter 8 of
		      Holger H. Hoos and Thomas Stuetzle, 
		      Stochastic Local Search---Foundations and Applications, 
		      Morgan Kaufmann Publishers, 2004.
		      or some of the papers online available from David S. Johnson.
*/
{
    long int c1, c2;             /* cities considered for an exchange */
    long int s_c1, s_c2;         /* successor cities of c1 and c2     */
    long int p_c1, p_c2;         /* predecessor cities of c1 and c2   */   
    long int pos_c1, pos_c2;     /* positions of cities c1, c2        */
    long int i, j, h, l;
    long int improvement_flag, improve_node, help, n_improves = 0, n_exchanges=0;
    long int h1=0, h2=0, h3=0, h4=0;
    long int radius;             /* radius of nn-search */
    long int gain = 0;
    long int *random_vector;
    long int *pos;               /* positions of cities in tour */ 
    long int *dlb;               /* vector containing don't look bits */ 
  
    pos = malloc(n * sizeof(long int));
    dlb = malloc(n * sizeof(long int));
    for ( i = 0 ; i < n ; i++ ) {
	pos[tour[i]] = i;
	dlb[i] = FALSE;
    }

    improvement_flag = TRUE;
    random_vector = generate_random_permutation( n );

    while ( improvement_flag ) {

	improvement_flag = FALSE;

	for (l = 0 ; l < n; l++) {

	    c1 = random_vector[l]; 
	    DEBUG ( assert ( c1 < n && c1 >= 0); )
		if ( dlb_flag && dlb[c1] )
		    continue;
	    improve_node = FALSE;
	    pos_c1 = pos[c1];
	    s_c1 = tour[pos_c1+1];
	    radius = instance.distance[c1][s_c1];

	    /* First search for c1's nearest neighbours, use successor of c1 */
	    for ( h = 0 ; h < nn_ls ; h++ ) {
		c2 = instance.nn_list[c1][h]; /* exchange partner, determine its position */
		if ( radius > instance.distance[c1][c2] ) {
		    s_c2 = tour[pos[c2]+1];
		    gain =  - radius + instance.distance[c1][c2] + 
			instance.distance[s_c1][s_c2] - instance.distance[c2][s_c2];
		    if ( gain < 0 ) {
			h1 = c1; h2 = s_c1; h3 = c2; h4 = s_c2; 
			improve_node = TRUE;
			goto exchange2opt;
		    }
		}
		else 
		    break;
	    }      
	    /* Search one for next c1's h-nearest neighbours, use predecessor c1 */
	    if (pos_c1 > 0)
		p_c1 = tour[pos_c1-1];
	    else 
		p_c1 = tour[n-1];
	    radius = instance.distance[p_c1][c1];
	    for ( h = 0 ; h < nn_ls ; h++ ) {
		c2 = instance.nn_list[c1][h];  /* exchange partner, determine its position */
		if ( radius > instance.distance[c1][c2] ) {
		    pos_c2 = pos[c2];
		    if (pos_c2 > 0)
			p_c2 = tour[pos_c2-1];
		    else 
			p_c2 = tour[n-1];
		    if ( p_c2 == c1 )
			continue;
		    if ( p_c1 == c2 )
			continue;
		    gain =  - radius + instance.distance[c1][c2] + 
			instance.distance[p_c1][p_c2] - instance.distance[p_c2][c2];
		    if ( gain < 0 ) {
			h1 = p_c1; h2 = c1; h3 = p_c2; h4 = c2; 
			improve_node = TRUE;
			goto exchange2opt;
		    }
		}
		else 
		    break;
	    }      
	    if (improve_node) {
	    exchange2opt:
		n_exchanges++;
		improvement_flag = TRUE;
		dlb[h1] = FALSE; dlb[h2] = FALSE;
		dlb[h3] = FALSE; dlb[h4] = FALSE;
		/* Now perform move */
		if ( pos[h3] < pos[h1] ) {
		    help = h1; h1 = h3; h3 = help;
		    help = h2; h2 = h4; h4 = help;
		}
		if ( pos[h3] - pos[h2] < n / 2 + 1) {
		    /* reverse inner part from pos[h2] to pos[h3] */
		    i = pos[h2]; j = pos[h3];
		    while (i < j) {
			c1 = tour[i];
			c2 = tour[j];
			tour[i] = c2;
			tour[j] = c1;
			pos[c1] = j;
			pos[c2] = i;
			i++; j--;
		    }
		}
		else {
		    /* reverse outer part from pos[h4] to pos[h1] */
		    i = pos[h1]; j = pos[h4];
		    if ( j > i )
			help = n - (j - i) + 1;
		    else 
			help = (i - j) + 1;
		    help = help / 2;
		    for ( h = 0 ; h < help ; h++ ) {
			c1 = tour[i];
			c2 = tour[j];
			tour[i] = c2;
			tour[j] = c1;
			pos[c1] = j;
			pos[c2] = i;
			i--; j++;
			if ( i < 0 )
			    i = n-1;
			if ( j >= n )
			    j = 0;
		    }
		    tour[n] = tour[0];
		}
	    } else {
		dlb[c1] = TRUE;
	    }
      
	}
	if ( improvement_flag ) {
	    n_improves++;
	}
    }
    free( random_vector );
    free( dlb );
    free( pos );
}

void two_h_opt_first( long int *tour ) 
/*    
      FUNCTION:       2-h-opt a tour
      INPUT:          pointer to the tour that undergoes local optimization
      OUTPUT:         none
      (SIDE)EFFECTS:  tour is 2-h-opt
      COMMENTS:       for details on 2-h-opt see J. L. Bentley. Fast algorithms for geometric 
                      traveling salesman problems. ORSA Journal on Computing, 
		      4(4):387--411, 1992. 
		      The neighbourhood is scanned in random order (this need 
                      not be the best possible choice). Concerning the speed-ups used 
		      here consult, for example, Chapter 8 of
		      Holger H. Hoos and Thomas Stuetzle, 
		      Stochastic Local Search---Foundations and Applications, 
		      Morgan Kaufmann Publishers, 2004.
		      or some of the papers online available from David S. Johnson.
*/
{

    long int c1, c2;         /* cities considered for an exchange */
    long int s_c1, s_c2;     /* successors of c1 and c2           */
    long int p_c1, p_c2;     /* predecessors of c1 and c2         */   
    long int pos_c1, pos_c2;     /* positions of cities c1, c2        */
    long int i, j, h, l;
    long int improvement_flag, improve_node;
    long int h1=0, h2=0, h3=0, h4=0, h5=0, help;
    long int radius;             /* radius of nn-search */
    long int gain = 0;
    long int *random_vector;
    long int two_move, node_move;

    long int *pos;               /* positions of cities in tour */ 
    long int *dlb;               /* vector containing don't look bits */ 
  
    pos = malloc(n * sizeof(long int));
    dlb = malloc(n * sizeof(long int));
    for ( i = 0 ; i < n ; i++ ) {
	pos[tour[i]] = i;
	dlb[i] = FALSE;
    }

    improvement_flag = TRUE;
    random_vector = generate_random_permutation( n );

    while ( improvement_flag ) {

	improvement_flag = FALSE; two_move = FALSE; node_move = FALSE;

	for (l = 0 ; l < n; l++) {

	    c1 = random_vector[l]; 
	    DEBUG ( assert ( c1 < n && c1 >= 0); )
		if ( dlb_flag && dlb[c1] )
		    continue;
	    improve_node = FALSE;
	    pos_c1 = pos[c1];
	    s_c1 = tour[pos_c1+1];
	    radius = instance.distance[c1][s_c1];

	    /* First search for c1's nearest neighbours, use successor of c1 */
	    for ( h = 0 ; h < nn_ls ; h++ ) {
		c2 = instance.nn_list[c1][h]; /* exchange partner, determine its position */
		if ( radius > instance.distance[c1][c2] ) {
		    pos_c2 = pos[c2];
		    s_c2 = tour[pos_c2+1];
		    gain =  - radius + instance.distance[c1][c2] + 
			instance.distance[s_c1][s_c2] - instance.distance[c2][s_c2];
		    if ( gain < 0 ) {
			h1 = c1; h2 = s_c1; h3 = c2; h4 = s_c2; 
			improve_node = TRUE; two_move = TRUE; node_move = FALSE;
			goto exchange;
		    }
		    if (pos_c2 > 0)
			p_c2 = tour[pos_c2-1];
		    else 
			p_c2 = tour[n-1];
		    gain = - radius + instance.distance[c1][c2] + instance.distance[c2][s_c1]
			+ instance.distance[p_c2][s_c2] - instance.distance[c2][s_c2]
			- instance.distance[p_c2][c2]; 
		    if ( c2 == s_c1 )
			gain = 0;
		    if ( p_c2 == s_c1 )
			gain = 0;

		    gain = 0;

		    if ( gain < 0 ) {
			h1 = c1; h2 = s_c1; h3 = c2; h4 = p_c2; h5 = s_c2; 
			improve_node = TRUE; node_move = TRUE; two_move = FALSE;
			goto exchange;
		    }
		}
		else 
		    break;
	    }      
	    /* Second search for c1's nearest neighbours, use predecessor c1 */
	    if (pos_c1 > 0)
		p_c1 = tour[pos_c1-1];
	    else 
		p_c1 = tour[n-1];
	    radius = instance.distance[p_c1][c1];
	    for ( h = 0 ; h < nn_ls ; h++ ) {
		c2 = instance.nn_list[c1][h];  /* exchange partner, determine its position */
		if ( radius > instance.distance[c1][c2] ) {
		    pos_c2 = pos[c2];
		    if (pos_c2 > 0)
			p_c2 = tour[pos_c2-1];
		    else 
			p_c2 = tour[n-1];
		    if ( p_c2 == c1 )
			continue;
		    if ( p_c1 == c2 )
			continue;
		    gain =  - radius + instance.distance[c1][c2] + 
			instance.distance[p_c1][p_c2] - instance.distance[p_c2][c2];
		    if ( gain < 0 ) {
			h1 = p_c1; h2 = c1; h3 = p_c2; h4 = c2; 
			improve_node = TRUE; two_move = TRUE; node_move = FALSE;
			goto exchange;
		    }
		    s_c2 = tour[pos[c2]+1];
		    gain = - radius + instance.distance[c2][c1] + instance.distance[p_c1][c2]
			+ instance.distance[p_c2][s_c2] - instance.distance[c2][s_c2]
			- instance.distance[p_c2][c2]; 
		    if ( p_c1 == c2 )
			gain = 0;
		    if ( p_c1 == s_c2 )
			gain = 0;

		    if ( gain < 0 ) {
			h1 = p_c1; h2 = c1; h3 = c2; h4 = p_c2; h5 = s_c2; 
			improve_node = TRUE; node_move = TRUE; two_move = FALSE;
			goto exchange;
		    }
		}
		else 
		    break;
	    }      
	exchange:
	    if (improve_node) {
		if ( two_move ) {
		    improvement_flag = TRUE;
		    dlb[h1] = FALSE; dlb[h2] = FALSE;
		    dlb[h3] = FALSE; dlb[h4] = FALSE;
		    /* Now perform move */
		    if ( pos[h3] < pos[h1] ) {
			help = h1; h1 = h3; h3 = help;
			help = h2; h2 = h4; h4 = help;
		    }
		    if ( pos[h3] - pos[h2] < n / 2 + 1) {
			/* reverse inner part from pos[h2] to pos[h3] */
			i = pos[h2]; j = pos[h3];
			while (i < j) {
			    c1 = tour[i];
			    c2 = tour[j];
			    tour[i] = c2;
			    tour[j] = c1;
			    pos[c1] = j;
			    pos[c2] = i;
			    i++; j--;
			}
		    }
		    else {
			/* reverse outer part from pos[h4] to pos[h1] */
			i = pos[h1]; j = pos[h4];
			if ( j > i )
			    help = n - (j - i) + 1;
			else 
			    help = (i - j) + 1;
			help = help / 2;
			for ( h = 0 ; h < help ; h++ ) {
			    c1 = tour[i];
			    c2 = tour[j];
			    tour[i] = c2;
			    tour[j] = c1;
			    pos[c1] = j;
			    pos[c2] = i;
			    i--; j++;
			    if ( i < 0 )
				i = n-1;
			    if ( j >= n )
				j = 0;
			}
			tour[n] = tour[0];
		    }
		} else if ( node_move ) {
		    improvement_flag = TRUE;
		    dlb[h1] = FALSE; dlb[h2] = FALSE; dlb[h3] = FALSE; 
		    dlb[h4] = FALSE; dlb[h5] = FALSE;
		    /* Now perform move */
		    if ( pos[h3] < pos[h1] ) {
			help = pos[h1] - pos[h3];
			i = pos[h3];
			for ( h = 0 ; h < help ; h++ ) {
			    c1 = tour[i+1];
			    tour[i] = c1;
			    pos[c1] = i;
			    i++;
			}
			tour[i] = h3;
			pos[h3] = i;
			tour[n] = tour[0];
		    } else {
			/* pos[h3] > pos[h1] */ 
			help = pos[h3] - pos[h1];
			/*  	    if ( help < n / 2 + 1) { */
			i = pos[h3];
			for ( h = 0 ; h < help - 1 ; h++ ) {
			    c1 = tour[i-1];
			    tour[i] = c1;
			    pos[c1] = i;
			    i--;
			}
			tour[i] = h3;
			pos[h3] = i;
			tour[n] = tour[0];
			/*  	  } */
		    }  
		} else {
		    fprintf(stderr,"this should never occur, 2-h-opt!!\n"); 
		    exit(0);
		}
		two_move = FALSE; node_move = FALSE;
	    } else {
		dlb[c1] = TRUE;
	    }
      
	}
    }
    free( random_vector );
    free( dlb );
    free( pos );
}

void three_opt_first( long int *tour )

/*    
      FUNCTION:       3-opt the tour
      INPUT:          pointer to the tour that is to optimize
      OUTPUT:         none
      (SIDE)EFFECTS:  tour is 3-opt
      COMMENT:        this is certainly not the best possible implementation of a 3-opt 
                      local search algorithm. In addition, it is very lengthy; the main 
		      reason herefore is that awkward way of making an exchange, where 
		      it is tried to copy only the shortest possible part of a tour. 
		      Whoever improves the code regarding speed or solution quality, please 
		      drop me the code at stuetzle no@spam informatik.tu-darmstadt.de
		      The neighbourhood is scanned in random order (this need 
                      not be the best possible choice). Concerning the speed-ups used 
		      here consult, for example, Chapter 8 of
		      Holger H. Hoos and Thomas Stuetzle, 
		      Stochastic Local Search---Foundations and Applications,