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