Revision 23229 trunk/extensions/extGraph/src/org/gvsig/graph/solvers/TspSolverAnnealing.java
TspSolverAnnealing.java | ||
---|---|---|
2 | 2 |
|
3 | 3 |
import java.util.Random; |
4 | 4 |
|
5 |
import org.gvsig.graph.core.GvFlag; |
|
6 |
|
|
5 | 7 |
import com.sun.org.apache.xpath.internal.operations.Mod; |
6 | 8 |
|
7 | 9 |
public class TspSolverAnnealing { |
10 |
static int T_INIT = 100; |
|
11 |
static double FINAL_T = 0.1; |
|
12 |
static double COOLING = 0.9; /* to lower down T (< 1) */ |
|
13 |
|
|
14 |
|
|
8 | 15 |
int n; |
9 | 16 |
int[] iorder; |
10 | 17 |
int[] jorder; |
11 | 18 |
float[] b = new float[4]; |
12 | 19 |
double[][] odMatrix; |
13 |
boolean g_bVolverOrigen;
|
|
20 |
boolean bReturnToOrigin;
|
|
14 | 21 |
int origenTSP, destinoTSP; |
15 | 22 |
int[] vTSP; |
16 | 23 |
|
17 | 24 |
|
18 | 25 |
|
19 |
static int T_INIT = 100; |
|
20 |
static double FINAL_T = 0.1; |
|
21 |
static double COOLING = 0.9; /* to lower down T (< 1) */ |
|
22 |
int TRIES_PER_T = 500*n; // TODO: PONER ESTO BIEN. |
|
23 |
int IMPROVED_PATH_PER_T = 60*n; // TODO: PONER ESTO BIEN. |
|
24 |
|
|
25 |
|
|
26 |
// static long A[56]= {-1}; |
|
27 |
// long *rand_fptr = A; |
|
28 |
|
|
29 |
// #define mod_diff(x,y) (((x)-(y))&0x7fffffff) |
|
30 |
// private long flipCycle() |
|
31 |
// { |
|
32 |
// long *ii,*jj; |
|
33 |
// for (ii = &A[1], jj = &A[32]; jj <= &A[55]; ii++, jj++) |
|
34 |
// *ii= mod_diff (*ii, *jj); |
|
35 |
// |
|
36 |
// for (jj = &A[1]; ii <= &A[55]; ii++, jj++) |
|
37 |
// *ii= mod_diff (*ii, *jj); |
|
38 |
// rand_fptr = &A[54]; |
|
39 |
// return A[55]; |
|
40 |
// } |
|
41 |
|
|
42 |
|
|
43 |
// typedef int Path[3]; /* specify how to change path */ |
|
44 |
|
|
45 | 26 |
/* |
46 | 27 |
* State vars |
47 | 28 |
*/ |
48 | 29 |
int verbose = 0; |
49 |
// Point *cities; |
|
50 |
// int *dist; |
|
30 |
private GvFlag[] flags = null; |
|
31 |
private boolean bVolverOrigen = false; |
|
32 |
private double distTotMin = Double.MAX_VALUE; |
|
51 | 33 |
|
52 |
// #define D(x,y) odMatrix[x][y] //dist[(x)*n+y] |
|
53 | 34 |
/* float calcula_dist_ordenacion(int v[], int bVolverOrigen) |
54 | 35 |
{ |
55 | 36 |
float dist, distTot; |
... | ... | |
58 | 39 |
distTot = odMatrix[origenTSP][v[0]]; // Origen al primer punto |
59 | 40 |
for (i = 0; i< numElemTSP-1;i++) |
60 | 41 |
{ |
61 |
//frmDocMap.Distancia v(i), v(i + 1), dist, tiempo |
|
62 | 42 |
dist = odMatrix[v[i]][v[i+1]]; |
63 | 43 |
distTot = distTot + dist; |
64 | 44 |
} |
65 | 45 |
|
66 | 46 |
// desde y hasta el almacen (distancia al primero y al ?ltimo |
67 |
// frmDocMap.Distancia idNodoOrigen, v(0), dist, tiempo |
|
68 | 47 |
if (bVolverOrigen) |
69 | 48 |
{ |
70 | 49 |
dist = odMatrix[v[numElemTSP-1]][origenTSP]; |
... | ... | |
84 | 63 |
int i; |
85 | 64 |
double len = 0; |
86 | 65 |
|
87 |
len = odMatrix[origenTSP][iorder[0]];
|
|
66 |
len = D(origenTSP, iorder[0]);
|
|
88 | 67 |
for (i = 0; i < n-1; i++) |
89 | 68 |
{ |
90 |
len += odMatrix[iorder[i]][iorder[i+1]];
|
|
69 |
len += D(iorder[i], iorder[i+1]);
|
|
91 | 70 |
} |
92 | 71 |
|
93 |
if (g_bVolverOrigen)
|
|
94 |
len += odMatrix[iorder[n-1]][origenTSP]; /* close path */
|
|
72 |
if (bReturnToOrigin)
|
|
73 |
len += D(iorder[n-1], origenTSP); /* close path */
|
|
95 | 74 |
else |
96 |
len += odMatrix[iorder[n-1]][destinoTSP]; /* close path */
|
|
75 |
len += D(iorder[n-1], destinoTSP); /* close path */
|
|
97 | 76 |
return (len); |
98 | 77 |
} |
99 | 78 |
|
100 |
|
|
101 |
/* |
|
102 |
* Prim's approximated TSP tour |
|
103 |
* See also [Cristophides'92] |
|
104 |
*/ |
|
105 |
// TODO: NO LO USO => Borrar esta funci?n. |
|
106 |
void findEulerianPath() |
|
107 |
{ |
|
108 |
int[] mst = new int[n]; |
|
109 |
int[] arc = new int[n]; |
|
110 |
int i, j = 0, k = 0, l, a; |
|
111 |
double d, maxd; |
|
112 |
|
|
113 |
double[] dis = new double[n]; |
|
114 |
|
|
115 |
maxd = Math.sqrt(b[1]-b[0])+ Math.sqrt(b[3]-b[2]); |
|
116 |
d = maxd; |
|
117 |
dis[0] = -1; |
|
118 |
for (i = 1; i < n; i++) |
|
119 |
{ |
|
120 |
dis[i] = odMatrix[i][0]; arc[i] = 0; |
|
121 |
if (d > dis[i]) |
|
122 |
{ |
|
123 |
d = dis[i]; |
|
124 |
j = i; |
|
125 |
} |
|
126 |
} |
|
127 |
|
|
128 |
/* |
|
129 |
* O(n^2) Minimum Spanning Trees by Prim and Jarnick |
|
130 |
* for graphs with adjacency matrix. |
|
131 |
*/ |
|
132 |
for (a = 0; a < n - 1; a++) |
|
133 |
{ |
|
134 |
mst[a] = j * n + arc[j]; /* join fragment j with MST */ |
|
135 |
dis[j] = -1; |
|
136 |
d = maxd; |
|
137 |
for (i = 0; i < n; i++) |
|
138 |
{ |
|
139 |
if (dis[i] >= 0) /* not connected yet */ |
|
140 |
{ |
|
141 |
if (dis[i] > odMatrix[i][j]) |
|
142 |
{ |
|
143 |
dis[i] = odMatrix[i][j]; |
|
144 |
arc[i] = j; |
|
145 |
} |
|
146 |
if (d > dis[i]) |
|
147 |
{ |
|
148 |
d = dis[i]; |
|
149 |
k = i; |
|
150 |
} |
|
151 |
} |
|
152 |
} |
|
153 |
j = k; |
|
154 |
} |
|
155 |
|
|
156 |
/* |
|
157 |
* Preorder Tour of MST |
|
158 |
*/ |
|
159 |
// #define VISITED(x) jorder[x] |
|
160 |
// #define NQ(x) arc[l++] = x |
|
161 |
// #define DQ() arc[--l] |
|
162 |
// #define EMPTY (l==0) |
|
163 |
|
|
164 |
for (i = 0; i < n; i++) jorder[i] = 0; |
|
165 |
k = 0; l = 0; d = 0; |
|
166 |
arc[l++] = 0; |
|
167 |
while (l != 0) |
|
168 |
{ |
|
169 |
i = arc[--l]; |
|
170 |
if (!(jorder[i] == 0)) |
|
171 |
{ |
|
172 |
iorder[k++] = i; |
|
173 |
jorder[i] = 1; |
|
174 |
for (j = 0; j < n - 1; j++) /* push all kids of i */ |
|
175 |
{ |
|
176 |
if (i == mst[j]%n) |
|
177 |
arc[l++] = mst[j]/n; |
|
178 |
} |
|
179 |
} |
|
180 |
} |
|
181 |
} |
|
182 | 79 |
|
183 | 80 |
int mod(int a, int b) { |
184 |
return (a % b); |
|
81 |
int aux = (a % b); |
|
82 |
if (aux < 0) |
|
83 |
aux = aux + b; |
|
84 |
return aux; |
|
185 | 85 |
} |
186 | 86 |
|
187 | 87 |
double D(int f, int t) { |
... | ... | |
206 | 106 |
|
207 | 107 |
// b va a ser el nuevo origen => sumamos su distancia a nuestro origen fijo |
208 | 108 |
double Dant, Dnuevo, Ddiff = 0; |
209 |
Dant = odMatrix[origenTSP][iorder[0]];
|
|
210 |
Dnuevo = odMatrix[origenTSP][b];
|
|
109 |
Dant = D(origenTSP, iorder[0]);
|
|
110 |
Dnuevo = D(origenTSP, b);
|
|
211 | 111 |
Ddiff = Dnuevo - Dant; |
212 | 112 |
|
213 | 113 |
// tambi?n hay que mirar la distancia al destino desde el ?ltimo punto |
214 | 114 |
int fin=n-1; |
215 |
if (g_bVolverOrigen)
|
|
115 |
if (bReturnToOrigin)
|
|
216 | 116 |
{ |
217 | 117 |
// p[2] va a ser el pr?ximo ?ltimo punto |
218 | 118 |
Dant = D(iorder[fin], origenTSP); |
... | ... | |
295 | 195 |
Ddiff = Dnuevo-Dant; |
296 | 196 |
} |
297 | 197 |
int fin = n-1; |
298 |
if (g_bVolverOrigen) // Miramos la distancia al cero
|
|
198 |
if (bReturnToOrigin) // Miramos la distancia al cero
|
|
299 | 199 |
{ |
300 | 200 |
// iorder[p[1]] o iorder[p[0]] va a acabar en la ?ltima posici?n |
301 | 201 |
if (p[0]==fin || p[1] == fin) // tambi?n hay que mirar la distancia al destino desde el ?ltimo punto |
... | ... | |
347 | 247 |
|
348 | 248 |
double annealing() |
349 | 249 |
{ |
350 |
int[] p; |
|
250 |
int[] p = new int[3];
|
|
351 | 251 |
int i=1, j, pathchg; |
352 | 252 |
int numOnPath, numNotOnPath; |
353 | 253 |
double pathlen, bestlen; |
354 | 254 |
double energyChange, T; |
355 |
|
|
255 |
|
|
356 | 256 |
/* |
357 | 257 |
* Set up first eulerian path iorder to be improved by |
358 | 258 |
* simulated annealing. |
... | ... | |
365 | 265 |
pathlen = pathLength(); // (iorder); |
366 | 266 |
bestlen = pathlen; |
367 | 267 |
Random rnd = new Random(); |
268 |
|
|
269 |
int TRIES_PER_T = 500*n; |
|
270 |
int IMPROVED_PATH_PER_T = 60*n; |
|
368 | 271 |
|
369 | 272 |
for (T = T_INIT; T > FINAL_T; T *= COOLING) /* annealing schedule */ |
370 | 273 |
{ |
... | ... | |
374 | 277 |
do { |
375 | 278 |
p[0] = rnd.nextInt(n); |
376 | 279 |
p[1] = rnd.nextInt(n); |
377 |
if (p[0] == p[1]) p[1] = mod(p[0]+1,n); /* non-empty path */ |
|
280 |
if (p[0] == p[1]) |
|
281 |
p[1] = mod(p[0]+1,n); /* non-empty path */ |
|
378 | 282 |
numOnPath = mod(p[1]-p[0],n) + 1; |
379 | 283 |
numNotOnPath = n - numOnPath; |
380 | 284 |
} while (numOnPath < 2 || numNotOnPath < 2) ; /* non-empty path */ |
381 | 285 |
|
382 |
|
|
286 |
double rreal = rnd.nextDouble(); |
|
383 | 287 |
if ((rnd.nextInt() % 2) == 0) /* threeWay */ |
384 | 288 |
{ |
385 | 289 |
do { |
386 | 290 |
p[2] = mod(rnd.nextInt(numNotOnPath)+p[1]+1,n); |
387 | 291 |
} while (p[0] == mod(p[2]+1,n)); /* avoids a non-change */ |
388 |
|
|
292 |
|
|
389 | 293 |
energyChange = getThreeWayCost (p); |
390 |
if ((energyChange < 0) || (RREAL < Math.exp(-energyChange/T)))
|
|
294 |
if ((energyChange < 0) || (rreal < Math.exp(-energyChange/T)))
|
|
391 | 295 |
{ |
392 | 296 |
pathchg++; |
393 | 297 |
pathlen += energyChange; |
... | ... | |
397 | 301 |
else /* path Reverse */ |
398 | 302 |
{ |
399 | 303 |
energyChange = getReverseCost (p); |
400 |
if ((energyChange < 0) || (RREAL < Math.exp(-energyChange/T)))
|
|
304 |
if ((energyChange < 0) || (rreal < Math.exp(-energyChange/T)))
|
|
401 | 305 |
{ |
402 | 306 |
pathchg++; |
403 | 307 |
pathlen += energyChange; |
... | ... | |
422 | 326 |
return bestlen; |
423 | 327 |
} |
424 | 328 |
|
329 |
|
|
330 |
public void setStops(GvFlag[] flags) { |
|
331 |
this.flags = flags; |
|
332 |
} |
|
333 |
|
|
334 |
|
|
335 |
/** |
|
336 |
* This function will use annealing only if numStops >= 6. If numStops < 6, it will |
|
337 |
* use direct calculation, ensuring the optimum result. |
|
338 |
* @return |
|
339 |
*/ |
|
340 |
public GvFlag[] calculate() { |
|
341 |
GvFlag[] orderedStops = new GvFlag[flags.length]; |
|
342 |
int numStops = flags.length; |
|
343 |
int numElemTSP; |
|
344 |
origenTSP = 0; |
|
345 |
destinoTSP = numStops-1; |
|
346 |
if (bVolverOrigen ) |
|
347 |
numElemTSP = numStops-1; |
|
348 |
else |
|
349 |
numElemTSP = numStops-2; |
|
350 |
|
|
351 |
n = numElemTSP; |
|
352 |
|
|
353 |
bReturnToOrigin = bVolverOrigen; |
|
354 |
|
|
355 |
iorder = new int[n]; |
|
356 |
jorder = new int[n]; |
|
357 |
|
|
358 |
|
|
359 |
|
|
360 |
vTSP = new int[numStops]; |
|
361 |
vTSP[0] = origenTSP; |
|
362 |
vTSP[numStops-1] = destinoTSP; |
|
363 |
if (numElemTSP > 0) |
|
364 |
{ |
|
365 |
// v = new int[numElemTSP]; |
|
366 |
for (int i=0; i< numElemTSP; i++) |
|
367 |
{ |
|
368 |
iorder[i] = i + 1; |
|
369 |
vTSP[i+1] = i + 1; // v[i]; |
|
370 |
} |
|
371 |
|
|
372 |
} |
|
373 |
|
|
374 |
/* identity permutation */ |
|
375 |
// for (i = 0; i < n; i++) iorder[i] = i+1; |
|
376 |
|
|
377 |
double distTotal=0.0; |
|
378 |
|
|
379 |
if (numElemTSP < 7) |
|
380 |
{ |
|
381 |
|
|
382 |
distTotMin = calculate_distance(iorder,bVolverOrigen); |
|
383 |
|
|
384 |
permut(iorder, numElemTSP, bVolverOrigen); |
|
385 |
|
|
386 |
distTotal = distTotMin ; |
|
387 |
} |
|
388 |
else |
|
389 |
{ |
|
390 |
distTotal = annealing(); |
|
391 |
} |
|
392 |
|
|
393 |
System.out.println("DistTotal = " + distTotal); |
|
394 |
System.out.print("new order = ["); |
|
395 |
for (int i=0; i < vTSP.length; i++) { |
|
396 |
orderedStops[i] = flags[vTSP[i]]; |
|
397 |
System.out.print(vTSP[i] + " "); |
|
398 |
} |
|
399 |
System.out.println("]"); |
|
400 |
return orderedStops; |
|
401 |
} |
|
402 |
|
|
403 |
|
|
404 |
double calculate_distance(int v[], boolean bVolverOrigen) |
|
405 |
{ |
|
406 |
double dist, distTot; |
|
407 |
int i; |
|
408 |
|
|
409 |
distTot = D(origenTSP, v[0]); // Origen al primer punto |
|
410 |
for (i = 0; i< n-1;i++) |
|
411 |
{ |
|
412 |
//frmDocMap.Distancia v(i), v(i + 1), dist, tiempo |
|
413 |
dist = D(v[i], v[i+1]); |
|
414 |
distTot = distTot + dist; |
|
415 |
} |
|
416 |
|
|
417 |
// desde y hasta el almacen (distancia al primero y al ?ltimo |
|
418 |
if (bVolverOrigen) |
|
419 |
{ |
|
420 |
dist = D(v[n-1], origenTSP); |
|
421 |
distTot = distTot + dist; |
|
422 |
} |
|
423 |
else |
|
424 |
{ |
|
425 |
dist = D(v[n-1], destinoTSP); |
|
426 |
distTot = distTot + dist; |
|
427 |
} |
|
428 |
return distTot; |
|
429 |
} |
|
430 |
|
|
431 |
void exchange(int v[], int p1, int p2) |
|
432 |
{ |
|
433 |
int aux; |
|
434 |
aux = v[p1]; |
|
435 |
v[p1] = v[p2]; |
|
436 |
v[p2] = aux; |
|
437 |
} |
|
438 |
|
|
439 |
void permut (int v[], int m, boolean bVolverOrigen) |
|
440 |
{ |
|
441 |
int i, j; |
|
442 |
double distTot; |
|
443 |
|
|
444 |
if (m > 0) |
|
445 |
for (i = 0; i < m; i++) |
|
446 |
{ |
|
447 |
exchange(v, i, m-1); |
|
448 |
permut (v, m-1, bVolverOrigen); |
|
449 |
exchange(v, i, m-1); |
|
450 |
} |
|
451 |
else |
|
452 |
{ |
|
453 |
// escribir_vector (v); |
|
454 |
// Aqu? tenemos el vector que buscamos |
|
455 |
// Calculamos su distancia y si es menor que la que tenemos guardada, guardamos el nuevo vector. |
|
456 |
distTot = calculate_distance(v, bVolverOrigen); |
|
457 |
|
|
458 |
if (distTot < distTotMin) |
|
459 |
{ |
|
460 |
distTotMin = distTot; |
|
461 |
for (j = 0; j< n; j++) |
|
462 |
vTSP[j+1] = v[j]; |
|
463 |
} |
|
464 |
// 'MuestraVector v |
|
465 |
|
|
466 |
} |
|
467 |
} |
|
468 |
|
|
469 |
|
|
470 |
public void setODMatrix(double[][] odMatrix) { |
|
471 |
this.odMatrix = odMatrix; |
|
472 |
} |
|
473 |
|
|
474 |
|
|
475 |
/** |
|
476 |
* @return true if the path is closed => will return to origin. |
|
477 |
* @see setReturnToOrigin |
|
478 |
*/ |
|
479 |
public boolean isClosed() { |
|
480 |
return bReturnToOrigin; |
|
481 |
} |
|
482 |
|
|
483 |
|
|
484 |
public void setReturnToOrigin(boolean returnToOrigin) { |
|
485 |
bReturnToOrigin = returnToOrigin; |
|
486 |
} |
|
487 |
|
|
425 | 488 |
} |
Also available in: Unified diff