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 }