1	/*
     2	 * Traveling sales person problem.
     3	 * example;
     4	 * > tsp 1.0 1.0 2.0 394
     5	 *
     6	 */
     7	
     8	#include <stdio.h>
     9	#include <math.h>
    10	#include <stdlib.h>
    11	
    12	#define CityNo    9
    13	
    14	/*#define rnd()     ( (double)rand() / 0x7fffffff )*/
    15	#define rnd()     ( (double)rand() / RAND_MAX )
    16	
    17	
    18	double distance[CityNo][CityNo];
    19	double unitout[CityNo][CityNo], Acoef, Bcoef, Dcoef;
    20	
    21	struct { double x, y; } cityxy[CityNo] =
    22	{ { 3.0, 4.0 }, { 2.0, 7.0 }, { 4.0, 7.0 }, { 5.0, 5.0 },
    23	    { 5.0, 3.0 }, { 4.0, 1.0 }, { 2.0, 1.0 }, { 1.0, 3.0 },
    24	    { 1.0, 5.0 } };
    25	
    26	int main( int argc, char *argv[] )
    27	{
    28	  int      no;
    29	  double   en;
    30	  unsigned seed;
    31	  void     initialize(), update_state(), display_state();
    32	  double   energy();
    33	
    34	  if ( argc <= 1 )
    35	    {
    36	      printf(" Usage: tsp Acoef Bcoef Dcoef seed\n");
    37	      exit(0);
    38	    }
    39	  Acoef = atof( argv[1] );
    40	  Bcoef = atof( argv[2] );
    41	  Dcoef = atof( argv[3] );
    42	  seed  = atoi( argv[4] );
    43	
    44	  printf("reading args: Acoef=%f, Bcoef=%f, Dcoef=%f, seed=%d\n",Acoef,Bcoef,Dcoef,seed);
    45	  srand( seed );
    46	
    47	  initialize();
    48	  display_state( 0 );
    49	  update_state();
    50	
    51	  for( no = 1; no < 50; no++ )
    52	    {
    53	      update_state();
    54	      display_state( no );
    55	      en = energy();
    56	      printf(" Energy = %f\n",en);
    57	    }
    58	
    59	  return( 0 );
    60	}
    61	
    62	void initialize()
    63	{
    64	  int    i, j;
    65	  double dtotal;
    66	
    67	  dtotal = 0.0;
    68	  for( i = 0; i < CityNo; i++ )
    69	    for( j = 0; j < CityNo; j++ )
    70	      {
    71		distance[i][j] = 0.1 *
    72		  sqrt( (cityxy[i].x - cityxy[j].x)*(cityxy[i].x - cityxy[j].x) +
    73		       (cityxy[i].y - cityxy[j].y)*(cityxy[i].y - cityxy[j].y));
    74		dtotal += distance[i][j];
    75	      }
    76	  
    77	  for( i = 0; i < CityNo; i++ )
    78	    for( j = 0; j < CityNo; j++)
    79	      distance[i][j] = 10.0*distance[i][j]/dtotal;
    80	
    81	  for( i = 0; i < CityNo; i++ )
    82	    for( j = 0; j < CityNo; j++ )
    83	      unitout[i][j] = rnd();
    84	}
    85	
    86	
    87	void update_state()
    88	{
    89	  int    i, j, n, m, jm, jp;
    90	  double un, unitin;
    91	  double aterm, bterm, dterm;
    92	
    93	  for( i = 0; i < CityNo; i++ )
    94	    for( j = 0; j < CityNo; j++ )
    95	      {
    96		aterm = bterm = dterm = 0.0;
    97		for( n = 0; n < CityNo; n++)
    98		  aterm += unitout[i][n];
    99		aterm = -Acoef*( aterm - unitout[i][j] );
   100	
   101		for( m = 0; m < CityNo; m++)
   102		  bterm += unitout[m][j];
   103		bterm = -Bcoef*( bterm - unitout[i][j] );
   104	
   105		if( j-1 == -1 ) jm = CityNo - 1;
   106		else jm = j - 1;
   107		if( j+1 == CityNo ) jp = 0;
   108		else jp = j + 1;
   109		for( m = 0; m < CityNo; m++)
   110		  dterm += distance[i][m]*(unitout[m][jp] + unitout[m][jm]);
   111		dterm = -Dcoef*dterm;
   112	
   113		unitin = aterm + bterm + dterm + Acoef + Bcoef;
   114		unitout[i][j] = 0.5*( 1.0 + tanh( unitin/0.5 ) );
   115	      }
   116	}
   117	
   118	void display_state( n )
   119	int n;
   120	{
   121	  int    i, j;
   122	
   123	  printf("    ###   Sequence     Cycle : %4d   ###\nCity ",n );
   124	  for( i = 0; i < CityNo; i++ )
   125	    printf(" %4d  ",i+1 );
   126	  printf("\n");
   127	  for( i = 0; i < CityNo; i++ )
   128	    {
   129	      printf("%4d  ",i+1 );
   130	      for( j = 0; j < CityNo; j++ )
   131		{
   132		  printf("%5.2f  ",unitout[i][j] );
   133		}
   134	      printf("\n");
   135	    }
   136	}
   137	
   138	double energy()
   139	{
   140	  int    i, j, m, jp, jm;
   141	  double term1, term2, term3;
   142	
   143	  term1 = term2 = term3 = 0;
   144	
   145	  for( i = 0; i < CityNo; i++ )
   146	    for( j = 0; j < CityNo; j++ )
   147	      for( m = 0; m < CityNo; m++ )
   148		if( m != j ) term1 += unitout[i][j]*unitout[i][m];
   149	  term1 = 0.5*Acoef*term1;
   150	
   151	  for( j = 0; j < CityNo; j++ )
   152	    for( i = 0; i < CityNo; i++ )
   153	      for( m = 0; m < CityNo; m++ )
   154		if( m != i ) term2 += unitout[i][j]*unitout[m][j];
   155	  term2 = 0.5*Bcoef*term2;
   156	
   157	  for( i = 0; i < CityNo; i++ )
   158	    for( j = 0; j < CityNo; j++ )
   159	      {
   160		if( j-1 == -1 ) jm = CityNo - 1;
   161		else jm = j - 1;
   162		if( j+1 == CityNo ) jp = 0;
   163		else jp = j + 1;
   164	       for( m = 0; m < CityNo; m++ )
   165		 term3 += distance[i][m]*unitout[i][j]*(unitout[m][jp]+unitout[m][jm]);
   166	      }
   167	  term3 = 0.5*Dcoef*term3;
   168	
   169	  return( term1 + term2 + term3 );
   170	}