#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <glpk.h>

typedef struct {
   int n;
   double *x;
   double *y;
   double **d;
} sites;

sites gensites ( int n )
{
   sites S;
   int i, j;

   S.n = n;
   S.x = (double *)malloc((n+1) * sizeof(double));
   S.y = (double *)malloc((n+1) * sizeof(double));
   for (i=1; i<=n; ++i) {
      S.x[i] = (double)rand() / (double)RAND_MAX;
      S.y[i] = (double)rand() / (double)RAND_MAX;
      printf("%18.16lf  %18.16lf\n", S.x[i], S.y[i]);
   }
   S.d = (double **)malloc((n+1) * sizeof(double *));
   S.d[0] = NULL;
   for (i=1; i<=n; ++i) S.d[i] = (double *)malloc((n+1) * sizeof(double));
   for (i=1; i<=n; ++i) {
      S.d[i][i] = 0;
      for (j=1; j<i; ++j) {
         S.d[i][j] = S.d[j][i] = sqrt((S.x[i]-S.x[j])*(S.x[i]-S.x[j]) + (S.y[i]-S.y[j])*(S.y[i]-S.y[j]));
      }
   }
   return S;
}

void tspsolve1 ( sites S )
{
   glp_prob *tsp;
   int n, i, j, k, l, I, J, nrows, ncols, ncoef;
   char name[32];
   int *ia, *ja;
   double *ar;
   int *next, *visited, *cycle, ncycles, clen, nattempt;

   n = S.n;

   next = (int *)malloc((n+1) * sizeof(int));
   visited = (int *)malloc((n+1) * sizeof(int));
   cycle = (int *)malloc((n*n+1) * sizeof(int));
   ncycles = nattempt = 0;

   while (1) {
      ++nattempt;
      printf("\n+++ Attempt No %d\n\n", nattempt);

      tsp = glp_create_prob();
      glp_set_prob_name(tsp, "TSP1");
      glp_set_obj_dir(tsp, GLP_MIN);

      ncols = n * n;
      glp_add_cols(tsp, ncols);
      for (i=1; i<=n; ++i) {
         for (j=1; j<=n; ++j) {
            sprintf(name, "x(%d,%d)", i, j);
            glp_set_col_name(tsp, n*(i-1)+j, name);
            glp_set_col_kind(tsp, n*(i-1)+j, GLP_BV);
            if (i == j) glp_set_col_bnds(tsp, n*(i-1)+j, GLP_FX, 0, 0);
            glp_set_obj_coef(tsp, n*(i-1)+j, S.d[i][j]);
         }
      }

      nrows = 2*n + n*(n-1)/2 + ncycles;
      glp_add_rows(tsp, nrows);
      for (i=1; i<=n; ++i) {
         sprintf(name, "from(%d)", i);
         glp_set_row_name(tsp, i, name);
         glp_set_row_bnds(tsp, i, GLP_FX, 1, 1);
         sprintf(name, "to(%d)", i);
         glp_set_row_name(tsp, n+i, name);
         glp_set_row_bnds(tsp, n+i, GLP_FX, 1, 1);
      }
      l = 2*n;
      for (i=1; i<=n; ++i) {
         for (j=1; j<i; ++j) {
            ++l;
            sprintf(name, "noloop(%d,%d)", i, j);
            glp_set_row_name(tsp, l, name);
            glp_set_row_bnds(tsp, l, GLP_DB, 0, 1);
         }
      }

      ncoef = nrows * ncols;
      ia = (int *)malloc((1+ncoef) * sizeof(int));
      ja = (int *)malloc((1+ncoef) * sizeof(int));
      ar = (double *)malloc((1+ncoef) * sizeof(double));
      l = 0;
      for (k=1; k<=n; ++k) {
         for (i=1; i<=n; ++i) {
            for (j=1; j<=n; ++j) {
               ++l;
               ia[l] = k;
               ja[l] = n*(i-1) + j;
               ar[l] = ((i != k) || (i == j)) ? 0 : 1;
               ++l;
               ia[l] = n+k;
               ja[l] = n*(i-1) + j;
               ar[l] = ((j != k) || (i == j)) ? 0 : 1;
            }
         }
      }
      k = 2*n;
      for (I=1; I<=n; ++I) {
         for (J=1; J<I; ++J) {
            ++k;
            for (i=1; i<=n; ++i) {
               for (j=1; j<=n; ++j) {
                  ++l;
                  ia[l] = k;
                  ja[l] = n*(i-1) + j;
                  if ((i == I) && (j == J)) ar[l] = 1;
                  else if ((i == J) && (j == I)) ar[l] = 1;
                  else ar[l] = 0;
               }
            }
         }
      }
      if (ncycles) {
         for (i=1; i<=n; ++i) visited[i] = 0;
         for (i=1; i<=n; ++i) {
            if (visited[i] == 0) {
               for (I=1; I<=n; ++I) for (J=1; J<=n; ++J) cycle[n*(I-1)+J] = 0;
               ++k;
               j = i;
               clen = 0;
               while (visited[j] == 0) {
                  visited[j] = 1;
                  cycle[n*(j-1)+next[j]] = 1;
                  j = next[j];
                  ++clen;
               }
               sprintf(name, "nocycle(%d)", i);
               glp_set_row_name(tsp, k, name);
               glp_set_row_bnds(tsp, k, GLP_FX, clen-1, clen-1);
               for (I=1; I<=n; ++I) {
                  for (J=1; J<=n; ++J) {
                     ++l;
                     ia[l] = k;
                     ja[l] = n*(I-1)+J;
                     ar[l] = cycle[n*(I-1)+J];
                  }
               }
            }
         }
      }

      glp_load_matrix(tsp, ncoef, ia, ja, ar);

      glp_simplex(tsp, NULL);

      printf("+++ Minimum simplex cost = %18.16lf\n", glp_get_obj_val(tsp));

      glp_intopt(tsp, NULL);

      printf("+++ Minimum mip cost     = %18.16lf\n\n", glp_mip_obj_val(tsp));
      for (i=1; i<=n; ++i) {
         for (j=1; j<=n; ++j) {
            k = glp_mip_col_val(tsp, n*(i-1)+j);
            if (k) next[i] = j;
         }
      }

      glp_delete_prob(tsp);
      free(ia); free(ja); free(ar);

      ncycles = 0;
      for (i=1; i<=n; ++i) visited[i] = 0;
      for (i=1; i<=n; ++i) {
         if (visited[i] == 0) {
            ++ncycles;
            printf("%-2d", i);
            j = i;
            while (visited[j] == 0) {
               visited[j] = 1;
               j = next[j];
               printf(" -> %-2d", j);
            }
            printf("\n");
         }
      }
      printf("\nSolution of Attempt %d has %d cycles\n", nattempt, ncycles);
      if (ncycles == 1) break;
   }

   free(next); free(visited); free(cycle);
}

void tspsolve2 ( sites S )
{
   glp_prob *tsp;
   int n, i, j, k, l, I, J, nrows, ncols, ncoef;
   char name[32];
   int *ia, *ja;
   double *ar;
   int *next, *visited, ncycles;

   n = S.n;

   next = (int *)malloc((n+1) * sizeof(int));
   visited = (int *)malloc((n+1) * sizeof(int));

   tsp = glp_create_prob();
   glp_set_prob_name(tsp, "TSP2");
   glp_set_obj_dir(tsp, GLP_MIN);

   ncols = n * n + n - 1;
   glp_add_cols(tsp, ncols);
   for (i=1; i<=n; ++i) {
      for (j=1; j<=n; ++j) {
         sprintf(name, "x(%d,%d)", i, j);
         glp_set_col_name(tsp, n*(i-1)+j, name);
         glp_set_col_kind(tsp, n*(i-1)+j, GLP_BV);
         if (i == j) glp_set_col_bnds(tsp, n*(i-1)+j, GLP_FX, 0, 0);
         glp_set_obj_coef(tsp, n*(i-1)+j, S.d[i][j]);
      }
   }
   for (i=2; i<=n; ++i) {
      sprintf(name, "u(%d)", i);
      glp_set_col_name(tsp, n*n+i-1, name);
      glp_set_col_kind(tsp, n*n+i-1, GLP_IV);
      glp_set_col_bnds(tsp, n*n+i-1, GLP_DB, 2, n);
      glp_set_obj_coef(tsp, n*n+i-1, 0);
   }

   nrows = 2*n + (n-1)*(n-2);
   glp_add_rows(tsp, nrows);
   for (i=1; i<=n; ++i) {
      sprintf(name, "from(%d)", i);
      glp_set_row_name(tsp, i, name);
      glp_set_row_bnds(tsp, i, GLP_FX, 1, 1);
      sprintf(name, "to(%d)", i);
      glp_set_row_name(tsp, n+i, name);
      glp_set_row_bnds(tsp, n+i, GLP_FX, 1, 1);
   }
   k = 2*n;
   for (i=2; i<=n; ++i) {
      for (j=2; j<=n; ++j) {
         if (i != j) {
            ++k;
            sprintf(name, "u(%d,%d)", i, j);
            glp_set_row_name(tsp, k, name);
            glp_set_row_bnds(tsp, k, GLP_UP, 0, n-1);
         }
      }
   }

   ncoef = nrows * ncols;
   ia = (int *)malloc((1+ncoef) * sizeof(int));
   ja = (int *)malloc((1+ncoef) * sizeof(int));
   ar = (double *)malloc((1+ncoef) * sizeof(double));
   l = 0;
   for (k=1; k<=n; ++k) {
      for (i=1; i<=n; ++i) {
         for (j=1; j<=n; ++j) {
            ++l;
            ia[l] = k;
            ja[l] = n*(i-1) + j;
            ar[l] = ((i != k) || (i == j)) ? 0 : 1;
            ++l;
            ia[l] = n+k;
            ja[l] = n*(i-1) + j;
            ar[l] = ((j != k) || (i == j)) ? 0 : 1;
         }
      }
      for (i=2; i<=n; ++i) {
         ++l;
         ia[l] = k;
         ja[l] = n*n + i - 1;
         ar[l] = 0;
         ++l;
         ia[l] = n+k;
         ja[l] = n*n + i - 1;
         ar[l] = 0;
      }
   }
   k = 2*n;
   for (i=2; i<=n; ++i) {
      for (j=2; j<=n; ++j) {
         if (i != j) {
            ++k;
            for (I=1; I<=n; ++I) {
               for (J=1; J<=n; ++J) {
                  ++l;
                  ia[l] = k;
                  ja[l] = n*(I-1)+J;
                  ar[l] = ((I == i) && (J == j)) ? n : 0;
               }
            }
            for (I=2; I<=n; ++I) {
               ++l;
               ia[l] = k;
               ja[l] = n*n + I - 1;
               if (I == i) ar[l] = 1;
               else if (I == j) ar[l] = -1;
               else ar[l] = 0;
            }
         }
      }
   }

   glp_load_matrix(tsp, ncoef, ia, ja, ar);

   glp_simplex(tsp, NULL);

   printf("+++ Minimum simplex cost = %18.16lf\n", glp_get_obj_val(tsp));

   glp_intopt(tsp, NULL);

   printf("+++ Minimum mip cost     = %18.16lf\n\n", glp_mip_obj_val(tsp));
   for (i=1; i<=n; ++i) {
      for (j=1; j<=n; ++j) {
         k = glp_mip_col_val(tsp, n*(i-1)+j);
         if (k) next[i] = j;
      }
   }

   glp_delete_prob(tsp);
   free(ia); free(ja); free(ar);

   ncycles = 0;
   for (i=1; i<=n; ++i) visited[i] = 0;
   for (i=1; i<=n; ++i) {
      if (visited[i] == 0) {
         ++ncycles;
         printf("%-2d", i);
         j = i;
         while (visited[j] == 0) {
            visited[j] = 1;
            j = next[j];
            printf(" -> %-2d", j);
         }
         printf("\n");
      }
   }
   printf("\nSolution has %d cycles\n", ncycles);

   free(next); free(visited);
}

int main ( int argc, char *argv[] )
{
   int n;
   sites S;

   srand((unsigned int)time(NULL));
   if (argc > 1) n = atoi(argv[1]); else scanf("%d", &n);
   printf("n = %d\n", n);
   S = gensites(n);

   glp_term_out(GLP_OFF);

   printf("\n++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n");
   printf("+++ Solving TSP by Method 1\n");
   printf("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n");
   tspsolve1(S);

   printf("\n++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n");
   printf("+++ Solving TSP by Method 2\n");
   printf("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\n");
   tspsolve2(S);

   exit(0);
}