// tsp.c  traveling salesman problem, NP Complete
//        test on 9 x,y random locations, 9! = 362,880
// arbitrary start and end, not closed loop
// will be different start and end for heuristic and complete

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "permute.h"

// problem data
static const int n = 9;
typedef struct {int id; double x; double y;} loc;
static loc locs[9];
static int soln[9];
static int avail[9];
static double total;
static int debug = 0;

static void used(int i, int na)
{ // delete i from avail that has n items by compressing
  int j, k;
  if(debug) printf("used(%d)\n", i);
  for(j=0; j<na; j++) if(avail[j]==i) {k=j; break;}
  for(j=k; j<na-1; j++) avail[j] = avail[j+1];
} // end used

static prsoln(int k, double dist)
{
  int i;

  if(debug) printf("debug k=%d, dist=%f, total=%f \n", k, dist, total);
  printf("soln  %d  ", soln[0]);
  for(i=1; i<k; i++) printf("%d  ", soln[i]);
  printf("\n");
  if(n-k>0)
  {
    printf("avail %d  ", avail[0]);
    for(i=1; i<n-k; i++) printf("%d  ", avail[i]);
    printf("\n");
  }
} // end prsoln

static int findca(int p, int na, double* dist)
{ // find closest available to p
  int i, j;
  double disti = 1.0E6;
  double xd, yd, dd;

  for(i=0; i<na; i++)
  { 
    if(debug) printf("compare p=%d to avail[%d]=%d locs[%d].id=%d\n", 
                     p, i, avail[i], avail[i], locs[avail[i]].id);
    xd = locs[p].x-locs[avail[i]].x;
    yd = locs[p].y-locs[avail[i]].y;
    dd = sqrt(xd*xd+yd*yd);
    if(dd<disti) { disti = dd; j = locs[avail[i]].id;} 
  }
  *dist = disti;
  return j;
} // end findca

int main(int argc, char * argv[])
{
  int i, j, k, j0, jk, last;
  long int ran1, ran2;
  double x, y, dist, dist0, distk, xd, yd, dd;
  double two31 = 1.0;
  int best[n]; // for NP Complete 9!

  printf("tsp.c running on %d locations\n", n);
  total = 0.0;
  for(i=0; i<31; i++) two31 *= 2.0;
  srandom(time(NULL));
  for(i=0; i<9; i++)
  {
    ran1 = random();
    x = (double)ran1/two31;
    ran2 = random();
    y = (double)ran2/two31;
    locs[i].id = i;
    locs[i].x = x;
    locs[i].y = y;
    printf("id=%d, x=%f, y=%f \n", locs[i].id, locs[i].x, locs[i].y);
    avail[i] = i; // not yet in soln
  }

  printf("heuristic solution running\n");
  total = 0.0;
  dist = 0.0;
  k = 0;
  // put in first point, find closest to first in soln
  soln[k] = 0;
  used(soln[k], n-k);
  k++;
  prsoln(k, dist);

  j = findca(soln[0], n-k, &dist);
  if(debug) printf("findca k=%d, j=%d, n-k=%d, dist=%f \n", k, j, n-k, dist); 
  soln[k] = j;
  used(soln[k], n-k);
  total += dist;
  k++;
  prsoln(k, dist);

  while(k<n)
  {
    j0 = findca(soln[0], n-k, &dist0);
    if(debug)
      printf("findca k=%d, j0=%d, n-k=%d, dist0=%f \n", k, j0, n-k, dist0); 
    jk = findca(soln[k-1], n-k, &distk);
    if(debug)
      printf("findca k=%d, jk=%d, n-k=%d, distk=%f \n", k, jk, n-k, distk); 
    if(distk<dist0)
    {
      soln[k] = jk;
      used(soln[k], n-k);
      dist = distk;
    }
    else
    { // move soln up, insert
      for(i=k; i>0; i--) soln[i] = soln[i-1]; 
      soln[0] = jk;
      used(soln[0], n-k);
      dist = dist0;
    }
    total += dist;
    k++;
    prsoln(k, dist);
  }
  printf("heuristic solution total distance=%f\n", total);

  // find minimum total distance for all permutations
  for(i=0; i<n; i++) soln[i] = i;
  avail[0] = n+1; // special initialization required for permuters
  last = 0;
  total = 1.0e6;

  while(!last)
  {
    // int permute(int n, int vector[], int permuters[]);
    last = permute(n, soln, avail);
    dist = 0.0;
    for(i=0; i<n-1; i++)
    {
      xd = locs[soln[i]].x-locs[soln[i+1]].x;
      yd = locs[soln[i]].y-locs[soln[i+1]].y;
      dd = sqrt(xd*xd+yd*yd);
      dist += dd;
    }
    if(dist<total)
    { 
      total = dist;
      for(i=0; i<n; i++) best[i] = soln[i]; 
    }
  }
  printf("best  %d  ", best[0]);
  for(i=1; i<k; i++) printf("%d  ", best[i]);
  printf("\n");
  printf("NP complete solution total distance=%f\n", total);

  printf("tsp.c finished\n");
  return 0;
} // end tsp.c