#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

/*
Approximately solve the traveling salesman problem (TSP) using
simulated annealing.
*/

/* The number of points. */
static int n = 100;

/* The coordinates of the points. */
static double* X;
static double* Y;

/* The starting temperature. */
static double T = 5;

/* The final temperature. */
static double TF = 0.05;


/* Calculate the distance travelled by visiting the points in X,Y in
   the order given in I. */
double func(const int* const I)
{
  double d = 0, u;
  int i;

  for (i=1; i<n; ++i)
    {
      u = (X[I[i]]-X[I[i-1]])*(X[I[i]]-X[I[i-1]]);
      u += (Y[I[i]]-Y[I[i-1]])*(Y[I[i]]-Y[I[i-1]]);
      d += sqrt(u);
    }

  return d;
}


/* Generate a step from I, put the result in Inew.  The step consists
   of a random number (rounded exponential with mean 2) of random
   swaps.
 */
void step(const int* const I,
	  int* const Inew)
{
  int r, i, j, k;

  memcpy(Inew, I, n*sizeof(int));

  for (r=0; r<1-2*log(drand48()); ++r)
    {
      /* Get a pair of distinct indices. */
      i = floor(n*drand48());
      while (1)
	{
	  j = floor(n*drand48());
	  if (i != j) break;
	}
      
      k = Inew[i];
      Inew[i] = Inew[j];
      Inew[j] = k;
    }
}


int main(int argc, char** argv)
{
  int i;
  int* I;
  int* Inew;
  int* IS;
  double sc, sc1;
  FILE* F;

  /* Storage for the points. */
  X = (double*)malloc(n*sizeof(double));
  Y = (double*)malloc(n*sizeof(double));

  /* Storage for the configurations. */
  I = (int*)malloc(n*sizeof(int));
  Inew = (int*)malloc(n*sizeof(int));

  /* Generate the points as iid uniform. */
  if (0)
    {
      for (i=0; i<n; ++i)
	{
	  X[i] = drand48();
	  Y[i] = drand48();
	}
    }

  /* Generate the points around the unit circle. */
  if (1)
    {
      for (i=0; i<n; ++i)
	{
	  X[i] = cos(2*M_PI*i/(double)n);
	  Y[i] = sin(2*M_PI*i/(double)n);
	}
    }


  /* Starting configuration. */
  for (i=0; i<n; ++i) I[i] = i;

  /* Score of starting configuration. */
  sc = func(I);

  while (T > TF)
    {
      step(I, Inew);
      sc1 = func(Inew);

      if (log(drand48()) < (sc-sc1)/T)
	{
	  /* Swap the pointers and update the score. */
	  IS = I;
	  I = Inew;
	  Inew = IS;
	  sc = sc1;

	  /* Cool a bit. */
	  T *= 0.9999;
	  printf("%f\t%f\n", T, sc);
	}

      break;
    }

  /* Output the final configuration. */
  F = fopen("anneal_out", "w");
  for (i=0; i<n; ++i)
    {
      fprintf(F, "%f\t%f\t%d\n", X[I[i]], Y[I[i]], I[i]);
    }
  fclose(F);

  return 0;
}