fork(3) download
copy 
  1. /*****************************************************************
  2.  Implementation of the robust taboo search of: E. Taillard
  3.  "Robust taboo search for the quadratic assignment problem", 
  4.  Parallel Computing 17, 1991, 443-455.
  5.  
  6.  Data file format: 
  7.   n, optimum solution value
  8.  (nxn) flow matrix,
  9.  (nxn) distance matrix
  10.  
  11.  Copyright : E. Taillard, 1990-2004
  12.  Standard C version with slight improvement regarding to
  13.  1991 version, E. Taillard, 14.03.2006
  14.  Compatibility: Unix and windows gcc, g++, bcc32.
  15.  This code can be freely used for non-commercial purpose.
  16.  Any use of this implementation or a modification of the code
  17.  must acknowledge the work of E. Taillard
  18.  
  19. ****************************************************************/
  20. #include <stdio.h>
  21. #include <stdlib.h>
  22.  
  23. const int infinite = 999999999;
  24. const int FALSE = 0;
  25. const int TRUE = 1;
  26.  
  27. typedef int*  type_vector;
  28. typedef int** type_matrix;
  29. int N[200]; /*solution*/
  30. int opt;
  31. double  somme_sol = 0.0;
  32. int to,count;
  33. long long int minimum_cost=2000000000;
  34.  
  35. /*************** L'Ecuyer random number generator ***************/
  36. double rando()
  37.  {
  38.   static int x10 = 12345, x11 = 67890, x12 = 13579, /* initial value*/
  39.              x20 = 24680, x21 = 98765, x22 = 43210; /* of seeds*/
  40.   const int m = 2147483647; const int m2 = 2145483479; 
  41.   const int a12= 63308; const int q12=33921; const int r12=12979; 
  42.   const int a13=-183326; const int q13=11714; const int r13=2883; 
  43.   const int a21= 86098; const int q21=24919; const int r21= 7417; 
  44.   const int a23=-539608; const int q23= 3976; const int r23=2071;
  45.   const double invm = 4.656612873077393e-10;
  46.   int h, p12, p13, p21, p23;
  47.   h = x10/q13; p13 = -a13*(x10-h*q13)-h*r13;
  48.   h = x11/q12; p12 = a12*(x11-h*q12)-h*r12;
  49.   if (p13 < 0) p13 = p13 + m; if (p12 < 0) p12 = p12 + m;
  50.   x10 = x11; x11 = x12; x12 = p12-p13; if (x12 < 0) x12 = x12 + m;
  51.   h = x20/q23; p23 = -a23*(x20-h*q23)-h*r23;
  52.   h = x22/q21; p21 = a21*(x22-h*q21)-h*r21;
  53.   if (p23 < 0) p23 = p23 + m2; if (p21 < 0) p21 = p21 + m2;
  54.   x20 = x21; x21 = x22; x22 = p21-p23; if(x22 < 0) x22 = x22 + m2;
  55.   if (x12 < x22) h = x12 - x22 + m; else h = x12 - x22;
  56.   if (h == 0) return(1.0); else return(h*invm);
  57.  }
  58.  
  59. /*********** return an integer between low and high *************/
  60. int unif(int low, int high)
  61.  {return low + (int)((double)(high - low + 1) * rando()) ;}
  62.  
  63. void transpose(int *a, int *b) {int temp = *a; *a = *b; *b = temp;}
  64.  
  65. int min(int a, int b) {if (a < b) return(a); else return(b);}
  66.  
  67. double cube(double x) {return x*x*x;}
  68.  
  69.  
  70. /*--------------------------------------------------------------*/
  71. /*       compute the cost difference if elements i and j        */
  72. /*         are transposed in permutation (solution) p           */
  73. /*--------------------------------------------------------------*/
  74. int compute_delta(int n, type_matrix a, type_matrix b,
  75.                    type_vector p, int i, int j)
  76.  {int d; int k;
  77.   d = (a[i][i]-a[j][j])*(b[p[j]][p[j]]-b[p[i]][p[i]]) +
  78.       (a[i][j]-a[j][i])*(b[p[j]][p[i]]-b[p[i]][p[j]]);
  79.   for (k = 0; k < n; k = k + 1) if (k!=i && k!=j)
  80.     d = d + (a[k][i]-a[k][j])*(b[p[k]][p[j]]-b[p[k]][p[i]]) +
  81.             (a[i][k]-a[j][k])*(b[p[j]][p[k]]-b[p[i]][p[k]]);
  82.   return(d);
  83.  }
  84.  
  85. /*--------------------------------------------------------------*/
  86. /*      Idem, but the value of delta[i][j] is supposed to       */
  87. /*    be known before the transposition of elements r and s     */
  88. /*--------------------------------------------------------------*/
  89. int compute_delta_part(type_matrix a, type_matrix b,
  90.                         type_vector p, type_matrix delta, 
  91.                         int i, int j, int r, int s)
  92.   {return(delta[i][j]+(a[r][i]-a[r][j]+a[s][j]-a[s][i])*
  93.      (b[p[s]][p[i]]-b[p[s]][p[j]]+b[p[r]][p[j]]-b[p[r]][p[i]])+
  94.      (a[i][r]-a[j][r]+a[j][s]-a[i][s])*
  95.      (b[p[i]][p[s]]-b[p[j]][p[s]]+b[p[j]][p[r]]-b[p[i]][p[r]]) );
  96.   }
  97.  
  98. void tabu_search(int n,                  /* problem size */
  99.                  type_matrix a,          /* flows matrix */
  100.                  type_matrix b,          /* distance matrix */
  101.                  type_vector best_sol,   /* best solution found */
  102.                  int *best_cost,         /* cost of best solution */
  103.                  int tabu_duration,      /* parameter 1 (< n^2/2) */
  104.                  int aspiration,         /* parameter 2 (> n^2/2)*/
  105.                  int nr_iterations)      /* number of iterations */ 
  106.  
  107.  
  108.  {type_vector p;                        /* current solution */
  109.   type_matrix delta;                    /* store move costs */
  110.   type_matrix tabu_list;                /* tabu status */
  111.   int current_iteration;                /* current iteration */
  112.   int current_cost;                     /* current sol. value */
  113.   int i, j, k, i_retained, j_retained;  /* indices */
  114.   int min_delta;                        /* retained move cost */
  115.   int autorized;                        /* move not tabu? */
  116.   int aspired;                          /* move forced? */
  117.   int already_aspired;                  /* in case many moves forced */
  118.  
  119.   /***************** dynamic memory allocation *******************/
  120.   p = (int*)calloc(n, sizeof(int));
  121.   delta = (int**)calloc(n,sizeof(int*));
  122.   for (i = 0; i < n; i = i+1) delta[i] = (int*)calloc(n, sizeof(int));
  123.   tabu_list = (int**)calloc(n,sizeof(int*));
  124.   for (i = 0; i < n; i = i+1) tabu_list[i] = (int*)calloc(n, sizeof(int));
  125.  
  126.  
  127.   /************** current solution initialization ****************/
  128.   for (i = 0; i < n; i = i + 1) p[i] = best_sol[i];
  129.  
  130.   /********** initialization of current solution value ***********/
  131.   /**************** and matrix of cost of moves  *****************/
  132.   current_cost = 0;
  133.   for (i = 0; i < n; i = i + 1) for (j = 0; j < n; j = j + 1)
  134.    {current_cost = current_cost + a[i][j] * b[p[i]][p[j]];
  135.     if (i < j) {delta[i][j] = compute_delta(n, a, b, p, i, j);};
  136.    };
  137.   *best_cost = current_cost;
  138.  
  139.   /****************** tabu list initialization *******************/
  140.   for (i = 0; i < n; i = i + 1) for (j = 0; j < n; j = j+1)
  141.     tabu_list[i][j] = -(n*i + j);
  142.  
  143.   /******************** main tabu search loop ********************/
  144.   for (current_iteration = 1; current_iteration <= nr_iterations;
  145.        current_iteration = current_iteration + 1)
  146.    {/** find best move (i_retained, j_retained) **/
  147.     i_retained = infinite;       /* in case all moves are tabu */
  148.     j_retained = infinite;
  149.     min_delta = infinite;
  150.     already_aspired = FALSE;
  151.  
  152.     for (i = 0; i < n-1; i = i + 1) 
  153.       for (j = i+1; j < n; j = j+1)
  154.        {autorized = (tabu_list[i][p[j]] < current_iteration) || 
  155.                     (tabu_list[j][p[i]] < current_iteration);
  156.  
  157.         aspired =
  158.          (tabu_list[i][p[j]] < current_iteration-aspiration)||
  159.          (tabu_list[j][p[i]] < current_iteration-aspiration)||
  160.          (current_cost + delta[i][j] < *best_cost);                
  161.  
  162.         if ((aspired && !already_aspired) || /* first move aspired*/
  163.            (aspired && already_aspired &&    /* many move aspired*/
  164.             (delta[i][j] < min_delta)   ) || /* => take best one*/
  165.            (!aspired && !already_aspired &&  /* no move aspired yet*/
  166.             (delta[i][j] < min_delta) && autorized))
  167.           {i_retained = i; j_retained = j;
  168.            min_delta = delta[i][j];
  169.            if (aspired) {already_aspired = TRUE;};
  170.           };
  171.        };
  172.  
  173.     if (i_retained == infinite) printf("All moves are tabu! \n"); 
  174.     else 
  175.      {/** transpose elements in pos. i_retained and j_retained **/
  176.       transpose(&p[i_retained], &p[j_retained]);
  177.       /* update solution value*/
  178.       current_cost = current_cost + delta[i_retained][j_retained];
  179.       /* forbid reverse move for a random number of iterations*/
  180.       tabu_list[i_retained][p[j_retained]] = 
  181.         current_iteration + (int)(cube(rando())*tabu_duration);
  182.       tabu_list[j_retained][p[i_retained]] = 
  183.          current_iteration + (int)(cube(rando())*tabu_duration);
  184.  
  185.       /* best solution improved ?*/
  186.       if (current_cost < *best_cost)
  187.        {*best_cost = current_cost;
  188.        for (k = 0; k < n; k = k+1) best_sol[k] = p[k];
  189.         printf("Solution of value: %d found at iter. %d\n", current_cost, current_iteration);
  190.        };
  191.            /* update matrix of the move costs*/
  192.       for (i = 0; i < n-1; i = i+1) for (j = i+1; j < n; j = j+1)
  193.         if (i != i_retained && i != j_retained && 
  194.             j != i_retained && j != j_retained)
  195.          {delta[i][j] = 
  196.             compute_delta_part(a, b, p, delta, 
  197.                                i, j, i_retained, j_retained);}
  198.         else
  199.          {delta[i][j] = compute_delta(n, a, b, p, i, j);};
  200.      };
  201.  
  202.    }; 
  203.   /* free memory*/
  204.   free(p);
  205.   for (i=0; i < n; i = i+1) free(delta[i]); free(delta);
  206.   for (i=0; i < n; i = i+1) free(tabu_list[i]); 
  207.   free(tabu_list);
  208. } /* tabu*/
  209.  
  210. void generate_random_solution(int n, type_vector  p)
  211.  {int i;
  212.   for (i = 0; i < n;   i++) p[i] = i;
  213.   for (i = 0; i < n-1; i++) transpose(&p[i], &p[unif(i, n-1)]);
  214.  }
  215.  
  216.  
  217.  
  218. void generate_solution(int n, type_vector p)
  219. {int i;
  220.  for(i = 0; i < n; i++) p[i] = N[i];
  221. } 
  222.  
  223. void define_initial_permutation(int n,type_vector p)
  224. {
  225.   int f, g;
  226.   int flag;
  227.   int used_num[200];
  228.  
  229.   /* 初期順列の生成 */
  230.  
  231.   for(f=0;f<n;f++)
  232.     {
  233.       used_num[f]=-1;
  234.     }
  235.  
  236.   for( f=0; f<n; f++)
  237.     {
  238.       do{
  239. 	p[f]=rand()%n;
  240. 	flag=0;
  241. 	for(g=0;g<n;g++)
  242. 	  {
  243. 	    if(used_num[g]==p[f])
  244. 	      flag=1;
  245. 	  }
  246.       }while(flag==1);
  247.       used_num[f]=p[f];
  248.     }
  249. }
  250.  
  251.  
  252.  
  253.  
  254.  
  255. int main()
  256.  {int n;                    /* problem size */
  257.   type_matrix a, b;         /* flows and distances matrices*/
  258.   type_vector solution;     /* solution (permutation) */
  259.   int cost;                 /* solution cost */
  260.   int pass_cost;
  261.    int nr_iterations,        /* number of tabu search iterations */
  262.       nr_resolutions, no_res; /* number of trials */
  263.  
  264.   FILE* data_file;       
  265.   char file_name[30];
  266.   int i, j;
  267.   char bidon[1000];
  268.  
  269.  
  270.   /************** read file name and problem size ***************/
  271.   printf("Data file name : \n");
  272.   scanf("%s",file_name); printf(file_name); printf("\n");
  273.   data_file = fopen(file_name,"r");
  274.   fscanf(data_file,"%d%d", &n, &opt);
  275.   fscanf(data_file,"%[^\n]", bidon);
  276.  
  277.   printf("Number of tabu search iterations:\n");
  278.   scanf("%d",&nr_iterations);
  279.   printf("Number of trials:\n");
  280.   scanf("%d",&nr_resolutions);
  281.  
  282.   /****************** dynamic memory allocation ******************/
  283.   solution = (int*)calloc(n, sizeof(int));
  284.  
  285.   a = (int**)calloc(n, sizeof(int*));
  286.   for (i = 0; i < n; i = i+1) a[i] = (int*)calloc(n, sizeof(int));
  287.   b = (int**)calloc(n,sizeof(int*));
  288.   for (i = 0; i < n; i = i+1) b[i] = (int*)calloc(n, sizeof(int));
  289.  
  290.   /************** read flows and distances matrices **************/
  291.   for (i = 0; i < n; i++) for (j = 0; j < n; j = j+1)
  292.     fscanf(data_file,"%d", &a[i][j]);
  293.  
  294.  
  295.   for (i = 0; i < n; i = i+1) for (j = 0; j < n; j = j+1)
  296.     fscanf(data_file,"%d", &b[i][j]);
  297.  
  298.   fclose(data_file);
  299.  
  300.   for(no_res = 1; no_res <= nr_resolutions; no_res++)
  301.     { printf("trials:%d\n",no_res);   
  302.  
  303.     if(no_res == 1){ 
  304.       srand(no_res);
  305.       define_initial_permutation(n,solution);
  306.       //  generate_random_solution(n, solution);
  307.     printf("first trial gone!\n");
  308.     } else 
  309.       {  generate_solution(n,solution);
  310.       printf("not first trial\n");
  311.       if(count == 10){
  312. 	printf("count is 10! \n");
  313.         count = 0;
  314.         srand(no_res);
  315. 	define_initial_permutation(n,solution);
  316. 	//      generate_random_solution(n,solution);
  317. }      }
  318.     tabu_search(n, a, b,                     /* problem data */
  319.                solution, &cost,              /* tabu search results */
  320.                8*n, n*n*5,                   /* parameters */
  321.                nr_iterations);               /* number of iterations */
  322.     printf("%d Solution found by tabu search:\n", cost);
  323.     if(cost == pass_cost)
  324.       count++;
  325.     else count = 0;
  326. 	  pass_cost = cost;
  327.  
  328.  
  329.      if(minimum_cost > cost)
  330.        minimum_cost = cost;
  331.      if (minimum_cost == opt)
  332.         to++;
  333.     for (i = 0; i < n; i = i+1){
  334.       N[i] = solution[i];
  335.     printf("%d ", solution[i]); 
  336.     }
  337.     printf("\n");
  338.     somme_sol += minimum_cost;
  339.    }
  340.    printf("trial:%d\n",nr_resolutions);
  341.    printf("Average cost: %f, average dev: %f\n", 
  342.           somme_sol/nr_resolutions, 100*(somme_sol/nr_resolutions - opt)/opt);
  343.    printf("Araival ave: %d, mimimum cost: %d\n" , to, minimum_cost);
  344.    printf(file_name);printf("\n");
  345.   free(solution);
  346.   for (i = 0; i < n; i = i+1) free(b[i]);
  347.   free(b);
  348.   for (i = 0; i < n; i = i+1) free(a[i]);
  349.   free(a);
  350.   fflush(stdin);
  351.   printf("Press return to end programme\n");
  352.   getchar();
  353.   return EXIT_SUCCESS;
  354.    }
Runtime error #stdin #stdout 0.01s 1860KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/FULhj
language:
C (gcc 8.3)
created:
12 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!