0votos

Held-Karp en C

por josejuan hace 5 años

Si sólo queremos el coste final, se puede reducir la memoria a O(n!/((n/2)!)^2) que es mucho menos que O(2^n) (aprox. un factor de 1/ln x). Aun se podría analizar el nº de MAXCOST que aparecen en un grafo total, xq no hace falta almacenarlos.

Implementar el algoritmo de Held-Karp para TSP simétrico.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
/* 
 
[josejuan@centella heldkarp]$ gcc -O2 -lm main.c 
[josejuan@centella heldkarp]$ time -f "%E, %M" ./a.out 25 
minCost := 23430.148438 
0:38.72, 509844 Kbytes 
 
*/ 
 
#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
 
#define MAXCOST 1000000.0f 
 
int comb(int n, int m) { 
    if(m) return (n * comb(n - 1, m - 1)) / m; 
    return 1; 
 
float heldKarpMinimumCost(int n, float (*cost)(int, int)) { 
 
    // max cost index 
    int maxComb = comb(n - 1, (n - 1) >> 1); 
 
    // las claves de combinación (sets de bits) ordenados de menor a mayor para búsqueda binaria 
    int *I1 = (int *) malloc(sizeof(int) * maxComb);                            // 2.704.156 * 4 
    int *I2 = (int *) malloc(sizeof(int) * maxComb);                            // 2.704.156 * 4 
 
    // costes del paso previo y actual para cada set y candidato de expansión 
    float *C1 = (float *) malloc(sizeof(float) * maxComb * (n - 1));            // 2.704.156 * 24 * 4 
    float *C2 = (float *) malloc(sizeof(float) * maxComb * (n - 1));            // 2.704.156 * 24 * 4 
                                                                                // ------------------------ 
                                                                                // 515 Mbytes 
 
 
    // I1 e incluso I2 se podrían eliminar (repitiendo cálculos) pero sólo reducirían como mucho 24M, 
    // de los 450M restantes únicamente se puede optimizar el hecho de que hay valores (¿cuantos?) que 
    // tienen coste máximo (no posibles) y usando un diccionario (en lugar de array) no consumirían 
    // memoria. 
 
    // número de combinaciones en cada registro 
    int c1; 
    int c2; 
 
    // estado inicial 
    c2 = 1; // sólo el nodo inicial 
    I2[0] = 0; // la clave es el set 0 (el primer nodo no cuenta) 
    // C2 (que será C1) no se usará la primera vez 
 
    int cmp(const void *u, const void *v) { 
        return *(int*)u - *(int*)v; 
    }; 
 
    int m; 
    for(m = 1; m < n; m++) { 
        { // intercambiamos registros 
            int   *I  = I1; float *C  = C1; int    c  = c1; 
                   I1 = I2;        C1 = C2;        c1 = c2; 
                   I2 = I ;        C2 = C ;        c2 = c ; 
        { // generamos combinaciones para el nuevo grupo 
            c2 = 0; 
 
            int set = (1 << m) - 1; 
            int limit = (1 << (n - 1)); 
            while (set < limit) { 
 
                I2[c2++] = set; 
 
                // Gosper's hack: 
                int c = set & -set; 
                int r = set + c; 
                set = (((r^set) >> 2) / c) | r; 
            int is, ic; 
            for(is = 0, ic = 0; is < c2; is++, ic += n - 1) { 
                int s = I2[is]; 
 
                int j = 1, mj = 1; 
                for(; j < n; j++, mj <<= 1) 
                    if(s & mj) { 
                        int s1 = s & ~mj; 
                        int *pis1 = (int*)bsearch(&s1, I1, c1, sizeof(int), cmp); // O(log (n!/(n/2)!(n/2)!)) 
                        int is1 = pis1 - I1; 
                        float minC = MAXCOST; // coste para k==1 
                        int k; 
                        for(k = 1; k < n; k++) 
                            if(k != j) { 
                                float ck = C1[is1 * (n - 1) + k - 1] + cost(k, j); 
                                if(ck < minC) 
                                    minC = ck; 
                        C2[ic + j - 1] = minC; 
                    } else { 
                        C2[ic + j - 1] = MAXCOST; 
 
    float minPathCost = MAXCOST; 
    { // en 2 el set total (con una sóla combinación) 
        int k; 
        for(k = 1; k < n; k++) { 
            float ck = C2[k - 1] + cost(k, 0); 
            if(ck < minPathCost) 
                minPathCost = ck; 
 
    free(I1); 
    free(I2); 
    free(C1); 
    free(C2); 
 
    return minPathCost; 
 
// ejemplo de uso de la función 
int main(int argc, char **argv) { 
 
 
    float sampleData [] = { 20833.3333f, 17100.0000f 
                          , 20900.0000f, 17066.6667f 
                          , 21300.0000f, 13016.6667f 
                          , 21600.0000f, 14150.0000f 
                          , 21600.0000f, 14966.6667f 
                          , 21600.0000f, 16500.0000f 
                          , 22183.3333f, 13133.3333f 
                          , 22583.3333f, 14300.0000f 
                          , 22683.3333f, 12716.6667f 
                          , 23616.6667f, 15866.6667f 
                          , 23700.0000f, 15933.3333f 
                          , 23883.3333f, 14533.3333f 
                          , 24166.6667f, 13250.0000f 
                          , 25149.1667f, 12365.8333f 
                          , 26133.3333f, 14500.0000f 
                          , 26150.0000f, 10550.0000f 
                          , 26283.3333f, 12766.6667f 
                          , 26433.3333f, 13433.3333f 
                          , 26550.0000f, 13850.0000f 
                          , 26733.3333f, 11683.3333f 
                          , 27026.1111f, 13051.9444f 
                          , 27096.1111f, 13415.8333f 
                          , 27153.6111f, 13203.3333f 
                          , 27166.6667f,  9833.3333f 
                          , 27233.3333f, 10450.0000f 
                          }; 
 
    int maxN = sizeof(sampleData) / (sizeof(float) << 1); 
 
    float *graph = (float *) malloc((sizeof(float) * maxN) << 5); // anchura fija de 5 bits (32 valores) 
 
    #define COST(a, b) graph[(a << 5) | b] 
    #define X(a) sampleData[a << 1] 
    #define Y(a) sampleData[(a << 1) + 1] 
    #define SQ(x) ((x)*(x)) 
 
    { // matriz de distancias 
        int a, b; 
        for(a = 0; a < maxN; a++) 
        for(b = 0; b < maxN; b++) 
            COST(a, b) = sqrt(SQ(X(a)-X(b))+SQ(Y(a)-Y(b))); 
 
    float cost(int a, int b) { 
        float c = COST(a, b); 
        return c; 
    }; 
 
    int N = atoi(argv[1]); 
 
    if(N < 0 || maxN < N) { 
        fprintf(stderr, "Argggggggggggggggggggggggggggg...!\n"); 
        return 1; 
 
    float minCost = heldKarpMinimumCost(N, cost); 
 
    printf("minCost := %f\n", minCost); 
 
    return 0; 
1 comentario
0votos

Escrito por josejuan hace 5 años

Una prueba empírica del nº de elementos con MAXCOST almacenados muestra que la mejora "sólo" es en un factor constante del 50% (es una serie alterna 0.5, 0.5+e1, 0.5, 0.5+e2, ... donde la serie e1, e2, ... tiende a cero).

[josejuan@centella heldkarp]$ for i in 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22; do ./a.out $i; done
LOWCOST	3	2	1	0.500000
LOWCOST	4	3	2	0.666667
LOWCOST	5	6	3	0.500000
LOWCOST	6	10	6	0.600000
LOWCOST	7	20	10	0.500000
LOWCOST	8	35	20	0.571429
LOWCOST	9	70	36	0.514286
LOWCOST	10	126	70	0.555556
LOWCOST	11	252	126	0.500000
LOWCOST	12	462	252	0.545455
LOWCOST	13	924	462	0.500000
LOWCOST	14	1716	924	0.538462
LOWCOST	15	3432	1717	0.500291
LOWCOST	16	6435	3432	0.533333
LOWCOST	17	12870	6435	0.500000
LOWCOST	18	24310	12870	0.529412
LOWCOST	19	48620	24311	0.500021
LOWCOST	20	92378	48620	0.526316
LOWCOST	21	184756	92379	0.500005
LOWCOST	22	352716	184756	0.523810
[josejuan@centella heldkarp]$ 

Comenta la solución

Tienes que identificarte para poder publicar tu comentario.