using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Text;
using System.IO;

namespace TPOpt
{
   class QAP
   {  /* Taboo search procedure for asymmetrical quadratic assignment problems */
      /* Author : Eric Taillard */
      /* Date : 94/11/23 */
      /* Implementation used in : E. Taillard, "Comparison of iterative searches */
      /* for the quadratic assignment problem", Publication 989, Center for */
      /* research on transportation, University of Montreal, Montreal, 1994 */
      
      const int infinite = 999999999;

      int n_max;
      int[,] a , b ;
      int[] solution, best_solution;
      double scout,cpu,cpu1,cpu2,cpusum,cpusum2,t2;
      int cost, i, n,init_sol_seed, meil_cout,
          no_res, nb_res, optimum, nb_iterations;    /* problem size */

      int best_cost, aspirating_iteration, nr_iteration_before_aspiration;
      int taboo_list_seed,current_taboo_list_size,lower_taboo_list_size, higher_taboo_list_size;
      int[,] taboo_list;

      DateTime t0;

      public QAP()
      {
      }

      public QAP (int n, int[,] d, int[,] f)
      {  this.n = n;
         a = d;
         b = f;
         optimum = 0;
         n_max   = n;
         solution      = new int[n];
         best_solution = new int[n];
         nb_iterations=n*2048;           /* di TS */
         nb_res=1;                       /* ripetizioni sullo stesso prob */
      }

      /// <summary>
      /// Distribuzione uniforme fra low e high
      /// </summary>
      /// <param name="seed"></param>
      /// <param name="low"></param>
      /// <param name="high"></param>
      /// <returns></returns>
      private int unif(int seed, int low, int high)
      {  int m = 2147483647, a = 16807, b = 127773, c = 2836;
         int kl;
         double value_0_1; // real double precision (coded on 64 bits)

         kl = seed / b;
         seed = a * (seed % b) - kl * c ;
         if (seed < 0) seed = seed + m ;
         value_0_1 = seed/(1.0*m);
         int res = low + (int) Math.Truncate((high - low + 1) * value_0_1);
         return res;
      }

      private void read_problem()
      {
         int i,j;
         StreamReader F = new StreamReader("NUG08.DAT");
         string line;
         string[] elem;
         char[] separators = new char[]{' '};

         optimum = Convert.ToInt32( F.ReadLine() );
         F.ReadLine();
         n_max = n = Convert.ToInt32( F.ReadLine() );

         a = new int[n_max, n_max];
         b = new int[n_max, n_max];
         solution = new int[n_max];
         best_solution = new int[n_max];

         F.ReadLine();
         nb_iterations=n*2048;           /* di TS */
         Trace.WriteLine(n+" "+optimum+" "+nb_res+" "+nb_iterations);
         nb_res=1;                       /* ripetizioni sullo stesso prob */

         for (i = 0; i<n; i++)
         {  line = F.ReadLine();
            elem = line.Split(separators,StringSplitOptions.RemoveEmptyEntries);
            for (j = 0; j<n; j++)
               a[i,j] = Convert.ToInt32( elem[j] );
         }
         F.ReadLine();
         for (i = 0; i<n; i++)
         {  line = F.ReadLine();
            elem = line.Split(separators,StringSplitOptions.RemoveEmptyEntries);
            for (j = 0; j<n; j++)
               b[i,j] = Convert.ToInt32( elem[j] );
         }

         F.Close();
      }

      private double clock()
      {
         TimeSpan tt = DateTime.Now - t0;
         return tt.TotalSeconds;
      }

      private void swap(ref int a, ref int b)
      {  int temp ;
         temp = a; 
         a = b; 
         b = temp;
      }

      private void improve_qap_solution(
              int n,  /* problem size */
              int[,] a, 
              int[,] b,   /* flows and distances matrices */
              int[] p,   /* initial solution -> best solution */
              ref int cost,  /* -> cost of the best solution */
              int optimum,
              int lower_taboo_list_size,                   
              int higher_taboo_list_size,                  
              int nr_iteration_before_aspiration,     /* parameters */
              int nr_iterations)
      {  int current_iteration, current_taboo_list_size, 
             aspirating_iteration, taboo_list_seed;

         taboo_list = new int[n_max,n_max];
         int[,] delta      = new int[n_max,n_max];    /* value of the difference between two solutions */
                  
         /* improve qap solution */
         initialize(taboo_list,p,delta);
         current_iteration = 0;
         do
         {  current_iteration = current_iteration+1;
            perform_one_move(current_iteration,delta,p,lower_taboo_list_size,higher_taboo_list_size);
            t2 = clock();
         } while ((best_cost > optimum) && (t2<15));
         cost = best_cost;
         p = best_solution;
      }

      /* computes the difference of solution values if units u and v are permuted */
      private void delta_full_computation(int i, int j, int[] p, int[,] delta)
      {  int k, sum;

         sum = 0;
         for (k = 0;k<n;k++)
         {  if (k!=i && k!=j) 
               sum = sum + a[k,i]*(b[p[k],p[j]]-b[p[k],p[i]]) +
                              a[k,j]*(b[p[k],p[i]]-b[p[k],p[j]]) + 
                              a[i,k]*(b[p[j],p[k]]-b[p[i],p[k]]) +
                              a[j,k]*(b[p[i],p[k]]-b[p[j],p[k]]);
               sum = sum + a[i,i]*(b[p[j],p[j]]-b[p[i],p[i]]) +
                           a[i,j]*(b[p[j],p[i]]-b[p[i],p[j]]) +
                           a[j,i]*(b[p[i],p[j]]-b[p[j],p[i]]) +
                           a[j,j]*(b[p[i],p[i]]-b[p[j],p[j]]);
            delta[i, j] = sum;
         }
      }

      /* idem but needs the value of delta[u,v] before the exchange of i and j  */
      private void delta_short_computation(int r, int s, int i, int j, int[,] delta, int[] p)
      {
            delta[i,j] = delta[i,j] + (a[r,i]-a[r,j]+a[s,j]-a[s,i])*
                  (b[p[s],p[i]]-b[p[s],p[j]]+b[p[r],p[j]]-b[p[r],p[i]]) +
                  (a[i,r]-a[j,r]+a[j,s]-a[i,s])*
                  (b[p[i],p[s]]-b[p[j],p[s]]+b[p[j],p[r]]-b[p[i],p[r]]);
      }

      private void initialize(int[,] taboo_list, int[] p, int[,] delta)
      {  int i, j;

         cost = 0;
         for (i=0;i<n;i++)
            for (j = 0;j<n;j++) 
            {
               taboo_list[i,j] = 0;
               cost = cost + a[i,j] * b[p[i],p[j]];
               if (i < j) 
               {
                  delta_full_computation(i, j, p, delta);
               }
            }
         Trace.WriteLine("Cost = "+cost);
         best_cost = cost;
         taboo_list_seed = 123456789;
         current_taboo_list_size = 
         unif(taboo_list_seed, lower_taboo_list_size, higher_taboo_list_size);

         nr_iteration_before_aspiration = 0;
         //lower_taboo_list_size = (int)Math.Round(0.9*n_max);
         //higher_taboo_list_size = (int)Math.Round(4.2*n_max);
      }
      
      private void find_best_move(ref int u, ref int v, int current_iteration, int[,] delta, int[] p)
      {  int i, j, delta_min;
         bool aspired, taboo;

         delta_min = infinite;
         u = infinite; /* in case no allowed moves exist */
         aspirating_iteration = current_iteration - nr_iteration_before_aspiration;
         for (i=0;i<n-1;i++) 
            for (j=i+1; j<n; j++) 
            {
               taboo = (taboo_list[i,p[j]] >= current_iteration) &&
                        (taboo_list[j,p[i]] >= current_iteration);
               aspired = ((taboo_list[i,p[j]] < aspirating_iteration) &&
                           (taboo_list[j,p[i]] < aspirating_iteration)) ||
                           ((cost + delta[i,j] < best_cost) && taboo);
               if (((delta[i,j] < delta_min) && !taboo) || aspired) 
               {
                  u = i; 
                  v = j;
                  if (aspired) 
                     delta_min = -infinite;
                  else 
                     delta_min = delta[i,j];
               }
            }
      }

      private void perform_one_move(int current_iteration, int[,] delta, int[] p, int lower_taboo_list_size, int higher_taboo_list_size)
      {
         int i, j, u=0, v=0;

         find_best_move(ref u,ref v,current_iteration,delta,p);
         if (u != infinite)
         {
            cost = cost + delta[u,v];
            taboo_list[u,p[u]] = current_iteration + current_taboo_list_size;
            taboo_list[v,p[v]] = current_iteration + current_taboo_list_size;
            swap(ref p[u], ref p[v]);
         }
         /* else trace.WriteLine('no allowed move !'); */
         if (cost < best_cost)
         {
            best_cost = cost;

            for(i=0;i<n;i++)
               best_solution[i] = p[i];
            cpu2 = clock();
            Trace.WriteLine(" iter "+current_iteration+" new cost "+best_cost+" t.cpu "+(cpu2 - cpu1));
         }

         for (i = 0; i<n-1;i++) 
            for (j = i+1;j<n;j++)
               if ((i!=u) && (i!=v) && (j!=u) && (j!=v))
                  delta_short_computation(u,v,i,j,delta,p);
               else 
                  delta_full_computation(i,j,p,delta);

         if (current_iteration % (2*higher_taboo_list_size) == 0)
            current_taboo_list_size =
               unif(taboo_list_seed, lower_taboo_list_size, higher_taboo_list_size);
   }

      public void qapTs ()
      {
         if(n<=0)
            read_problem();
         init_sol_seed = 123456789;
         scout = 0; 
         meil_cout = infinite;
         t0 = DateTime.Now;
         cpusum=0;
         cpusum2=0;
         for (no_res=0; no_res<nb_res;no_res++) 
         {
            Trace.WriteLine("=====> Repetition n."+no_res);
            cpu1=clock();
            Trace.WriteLine("cpu1="+cpu1);
            for (i = 0;i<n;i++)
               solution[i] = i;
            for (i = 0;i<n-1;i++)
               swap(ref solution[i], ref solution[unif(init_sol_seed, i, n)]);

            improve_qap_solution(n, a, b, solution, ref cost, optimum,
                              (int) Math.Truncate(0.9*n), (int) Math.Truncate(4.1*n) + 4, 3*n*n, nb_iterations);
            cpu2=clock();
            Trace.WriteLine("cpu2="+cpu2);
            cpu=cpu2-cpu1;
            if(cpu>0)  
               Trace.WriteLine("(1) cpu="+cpu+" cpu2="+cpu2+" cpu1="+cpu1+" cost:"+cost);
            cpusum=cpusum+cpu;
            scout = scout + cost;
            if(cost<meil_cout) 
               meil_cout=cost;
         }
         int cst = 0;
         for(i=0;i<n;i++)
         {  Trace.Write(best_solution[i]+" ");
            for(int j=0;j<n;j++)
               cst += a[i,j]*b[best_solution[i],best_solution[j]];
         }
         Trace.WriteLine("costo "+cst);
         Trace.WriteLine("Results: avg "+scout/nb_res+" best "+meil_cout+" t.avg (sec) "+cpusum/(1.00*nb_res));
      }
   }
}