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

#define NMAX 10000
#define True 1
#define False 0


double length(int N, double *x, double *y);
float ran2(long *idum);

main(int argc,char **argv)
{
  int i,N,i1,i2,n,accept,niter,naccept,nreject,circle;
  long seed=113097;
  
  double T,T0,dT,dT0,l1,l2,P,lglobalmin,lorig,changedT1,changedT2,phi;
  double x[NMAX],y[NMAX],xtmp,ytmp;

  char buf[120],city[NMAX][80],citytmp[80];
  FILE *fp;

  if (argc<3) {
    printf("Usage: %s Npoints T0 dT [changedT1 [changedT2]] \n",argv[0]);
    printf("If Npoints=0 read in coordinates from stdin \n");
    printf("If Npoints<0 generate coordinates on a circle\n");
    exit(0);
  }

  sscanf(argv[1],"%d",&N);
  sscanf(argv[2],"%lg",&T0);
  sscanf(argv[3],"%lg",&dT);

  circle=False;
  if (N<0) {
    N=-N;
    circle=True;
  }

  changedT1=-1.0;
  changedT2=-1.0;
  if (argc>4) sscanf(argv[4],"%lg",&changedT1);
  if (argc>5) sscanf(argv[5],"%lg",&changedT2);
  printf("changedT %g %g\n",changedT1,changedT2);
  
  if (N>NMAX) { 
    printf("Error: N too large, increase NMAX\n"); exit(0);
  }

  for(i=0;i<NMAX;i++) { city[i][0]='\0'; }

  if (N==0) {
    i=0;
    while(fgets(buf,120,stdin) != NULL) {
      n=sscanf(buf,"%lg %lg %79s",&(x[i]),&(y[i]),city[i]);
      if (n>=2) i++;
    }
    N=i;
    printf("Read in %d cities\n",N);
  }
  else {
    for(i=0;i<N;i++) {
      if (circle) {
	phi=2*3.1415927*ran2(&seed);
	x[i]=cos(phi);
	y[i]=sin(phi);
      }
      else {
	x[i]=ran2(&seed);
	y[i]=ran2(&seed);
      }
    }
  }

  fp=fopen("traveling_origpath","w");
  for (i=0;i<N;i++) {
    fprintf(fp,"%g %g\n",x[i],y[i]);
  }
  fprintf(fp,"%g %g %s\n",x[0],y[0],city[0]);
  fclose(fp);
  
  T=T0;
  dT0=dT;
  printf("Solving the traveler problem with T0=%g dT=%g\n",T0,dT);
  l1=length(N,x,y);
  lorig=l1;
  printf("Initial length %g\n",l1);    
  
  lglobalmin=l1;
  naccept=0; nreject=0;
  while (T > 0.0) {

    /* Swap 2 positions */
    i1=(int) (ran2(&seed)*N); if (i1>N-1) i1=N;
    i2=(int) (ran2(&seed)*N); if (i2>N-1) i2=N;
    
    xtmp=x[i1];  ytmp=y[i1];  strncpy(citytmp,city[i1],80);
    x[i1]=x[i2]; y[i1]=y[i2]; strncpy(city[i1],city[i2],80);
    x[i2]=xtmp;  y[i2]=ytmp;  strncpy(city[i2],citytmp,80);

    l2=length(N,x,y);

    /* Metropolis algoritm */
    accept=False;
    if (l2 < l1) accept=True;
    else {
      P=exp(-(l2-l1)/T);
      if (ran2(&seed) < P) accept=True;
    }

    if (accept) {
      l1=l2;
      naccept++;
      if (naccept%1000==0) printf("New length %10g global min %10g at T=%10g dt=%g\n",l1,lglobalmin,T,dT);
    }
    else {
      nreject++;
      /* Swap back */
      xtmp=x[i1];  ytmp=y[i1];   strncpy(citytmp,city[i1],80);
      x[i1]=x[i2]; y[i1]=y[i2];  strncpy(city[i1],city[i2],80);
      x[i2]=xtmp;  y[i2]=ytmp;   strncpy(city[i2],citytmp,80);
    }

    if (l1 < lglobalmin) lglobalmin=l1;

    /* Simulated annealing: lower T continuously */
    T=T-dT;
    if (T<changedT1) dT=dT0/10.0;
    if (T<changedT2) dT=dT0/100.0;

    niter++;

  }

  printf("New length %10g global min %10g at T=%10g dT=%g\n",l1,lglobalmin,T,dT);  

  fp=fopen("traveling_path","w");
  for (i=0;i<N;i++) {
    fprintf(fp,"%g %g %s\n",x[i],y[i],city[i]);
  }
  fprintf(fp,"%g %g %s\n",x[0],y[0],city[0]);
  fclose(fp);

  printf("Result: original length %g final %g naccept %d nreject %d\n",lorig,l1,naccept,nreject);

}

double length(int N, double *x, double *y)
{
  int i;
  double dl,l;

  l=0.0;
  for (i=1;i<N;i++) {
    dl=sqrt( (x[i]-x[i-1])*(x[i]-x[i-1]) + (y[i]-y[i-1])*(y[i]-y[i-1]));
    l+=dl;
  }
  dl=sqrt( (x[N-1]-x[0])*(x[N-1]-x[0]) + (y[N-1]-y[0])*(y[N-1]-y[0]));
  l+=dl;
  

  return l;

}




#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

float ran2(long *idum)
{
	int j;
	long k;
	static long idum2=123456789;
	static long iy=0;
	static long iv[NTAB];
	float temp;

	if (*idum <= 0) {
		if (-(*idum) < 1) *idum=1;
		else *idum = -(*idum);
		idum2=(*idum);
		for (j=NTAB+7;j>=0;j--) {
			k=(*idum)/IQ1;
			*idum=IA1*(*idum-k*IQ1)-k*IR1;
			if (*idum < 0) *idum += IM1;
			if (j < NTAB) iv[j] = *idum;
		}
		iy=iv[0];
	}
	k=(*idum)/IQ1;
	*idum=IA1*(*idum-k*IQ1)-k*IR1;
	if (*idum < 0) *idum += IM1;
	k=idum2/IQ2;
	idum2=IA2*(idum2-k*IQ2)-k*IR2;
	if (idum2 < 0) idum2 += IM2;
	j=iy/NDIV;
	iy=iv[j]-idum2;
	iv[j] = *idum;
	if (iy < 1) iy += IMM1;
	if ((temp=AM*iy) > RNMX) return RNMX;
	else return temp;
}
#undef IM1
#undef IM2
#undef AM
#undef IMM1
#undef IA1
#undef IA2
#undef IQ1
#undef IQ2
#undef IR1
#undef IR2
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX