Dump pbqp solution to file...
[libfirm] / heuristical.c
1 #include "adt/array.h"
2 #include "assert.h"
3 #include "error.h"
4
5 #include "bucket.h"
6 #include "heuristical.h"
7 #include "html_dumper.h"
8 #include "kaps.h"
9 #include "matrix.h"
10 #include "pbqp_edge.h"
11 #include "pbqp_edge_t.h"
12 #include "pbqp_node.h"
13 #include "pbqp_node_t.h"
14 #include "vector.h"
15
16 static pbqp_edge **edge_bucket;
17 static pbqp_node **node_buckets[4];
18 static pbqp_node **reduced_bucket = NULL;
19 static int         buckets_filled = 0;
20
21 static int dump = 0;
22
23 /* Forward declarations. */
24 static void apply_Brute_Force(pbqp *pbqp);
25
26 static void insert_into_edge_bucket(pbqp_edge *edge)
27 {
28         if (edge_bucket_contains(edge_bucket, edge)) {
29                 /* Edge is already inserted. */
30                 return;
31         }
32
33         edge_bucket_insert(&edge_bucket, edge);
34 }
35
36 static void init_buckets(void)
37 {
38         int i;
39
40         edge_bucket_init(&edge_bucket);
41         node_bucket_init(&reduced_bucket);
42
43         for (i = 0; i < 4; ++i) {
44                 node_bucket_init(&node_buckets[i]);
45         }
46 }
47
48 static void free_buckets(void)
49 {
50         int i;
51
52         for (i = 0; i < 4; ++i) {
53                 node_bucket_free(&node_buckets[i]);
54         }
55
56         edge_bucket_free(&edge_bucket);
57         node_bucket_free(&reduced_bucket);
58
59         buckets_filled = 0;
60 }
61
62 static void fill_node_buckets(pbqp *pbqp)
63 {
64         unsigned node_index;
65         unsigned node_len;
66
67         assert(pbqp);
68         node_len = pbqp->num_nodes;
69
70         for (node_index = 0; node_index < node_len; ++node_index) {
71                 unsigned   degree;
72                 pbqp_node *node = get_node(pbqp, node_index);
73
74                 if (!node) continue;
75
76                 degree = pbqp_node_get_degree(node);
77
78                 /* We have only one bucket for nodes with arity >= 3. */
79                 if (degree > 3) {
80                         degree = 3;
81                 }
82
83                 node_bucket_insert(&node_buckets[degree], node);
84         }
85
86         buckets_filled = 1;
87 }
88
89 static void normalize_towards_source(pbqp *pbqp, pbqp_edge *edge)
90 {
91         pbqp_matrix    *mat;
92         pbqp_node      *src_node;
93         pbqp_node      *tgt_node;
94         vector         *src_vec;
95         vector         *tgt_vec;
96         int             src_len;
97         int             tgt_len;
98         int             src_index;
99
100         assert(pbqp);
101         assert(edge);
102
103         src_node = edge->src;
104         tgt_node = edge->tgt;
105         assert(src_node);
106         assert(tgt_node);
107
108         src_vec = src_node->costs;
109         tgt_vec = tgt_node->costs;
110         assert(src_vec);
111         assert(tgt_vec);
112
113         src_len = src_vec->len;
114         tgt_len = tgt_vec->len;
115         assert(src_len > 0);
116         assert(tgt_len > 0);
117
118         mat = edge->costs;
119         assert(mat);
120
121         /* Normalize towards source node. */
122         for (src_index = 0; src_index < src_len; ++src_index) {
123                 num min = pbqp_matrix_get_row_min(mat, src_index, tgt_vec);
124
125                 if (min != 0) {
126                         if (src_vec->entries[src_index].data == INF_COSTS) {
127                                 pbqp_matrix_set_row_value(mat, src_index, 0);
128                         } else {
129                                 pbqp_matrix_sub_row_value(mat, src_index, tgt_vec, min);
130                         }
131                         src_vec->entries[src_index].data = pbqp_add(
132                                         src_vec->entries[src_index].data, min);
133
134                         if (min == INF_COSTS) {
135                                 unsigned edge_index;
136                                 unsigned edge_len = pbqp_node_get_degree(src_node);
137
138                                 for (edge_index = 0; edge_index < edge_len; ++edge_index) {
139                                         pbqp_edge *edge_candidate = src_node->edges[edge_index];
140                                         if (edge_candidate != edge) {
141                                                 insert_into_edge_bucket(edge_candidate);
142                                         }
143                                 }
144                         }
145                 }
146         }
147 }
148
149 static void normalize_towards_target(pbqp *pbqp, pbqp_edge *edge)
150 {
151         pbqp_matrix    *mat;
152         pbqp_node      *src_node;
153         pbqp_node      *tgt_node;
154         vector         *src_vec;
155         vector         *tgt_vec;
156         int             src_len;
157         int             tgt_len;
158         int             tgt_index;
159
160         assert(pbqp);
161         assert(edge);
162
163         src_node = edge->src;
164         tgt_node = edge->tgt;
165         assert(src_node);
166         assert(tgt_node);
167
168         src_vec = src_node->costs;
169         tgt_vec = tgt_node->costs;
170         assert(src_vec);
171         assert(tgt_vec);
172
173         src_len = src_vec->len;
174         tgt_len = tgt_vec->len;
175         assert(src_len > 0);
176         assert(tgt_len > 0);
177
178         mat = edge->costs;
179         assert(mat);
180
181         for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
182                 num min = pbqp_matrix_get_col_min(mat, tgt_index, src_vec);
183
184                 if (min != 0) {
185                         if (tgt_vec->entries[tgt_index].data == INF_COSTS) {
186                                 pbqp_matrix_set_col_value(mat, tgt_index, 0);
187                         } else {
188                                 pbqp_matrix_sub_col_value(mat, tgt_index, src_vec, min);
189                         }
190                         tgt_vec->entries[tgt_index].data = pbqp_add(
191                                         tgt_vec->entries[tgt_index].data, min);
192
193                         if (min == INF_COSTS) {
194                                 unsigned edge_index;
195                                 unsigned edge_len = pbqp_node_get_degree(tgt_node);
196
197                                 for (edge_index = 0; edge_index < edge_len; ++edge_index) {
198                                         pbqp_edge *edge_candidate = tgt_node->edges[edge_index];
199                                         if (edge_candidate != edge) {
200                                                 insert_into_edge_bucket(edge_candidate);
201                                         }
202                                 }
203                         }
204                 }
205         }
206 }
207
208 static void reorder_node(pbqp_node *node)
209 {
210         unsigned    degree     = pbqp_node_get_degree(node);
211         /* Assume node lost one incident edge. */
212         unsigned    old_degree = degree + 1;
213
214         if (!buckets_filled) return;
215
216         /* Same bucket as before */
217         if (degree > 2) return;
218
219         if (!node_bucket_contains(node_buckets[old_degree], node)) {
220                 /* Old arity is new arity, so we have nothing to do. */
221                 assert(node_bucket_contains(node_buckets[degree], node));
222                 return;
223         }
224
225         /* Delete node from old bucket... */
226         node_bucket_remove(&node_buckets[old_degree], node);
227
228         /* ..and add to new one. */
229         node_bucket_insert(&node_buckets[degree], node);
230 }
231
232 static void check_melting_possibility(pbqp *pbqp, pbqp_edge *edge)
233 {
234         pbqp_matrix    *mat;
235         pbqp_node      *src_node;
236         pbqp_node      *tgt_node;
237         vector         *src_vec;
238         vector         *tgt_vec;
239         int             src_len;
240         int             tgt_len;
241         int             src_index;
242         int             tgt_index;
243
244         assert(pbqp);
245         assert(edge);
246
247         src_node = edge->src;
248         tgt_node = edge->tgt;
249         assert(src_node);
250         assert(tgt_node);
251
252         src_vec = src_node->costs;
253         tgt_vec = tgt_node->costs;
254         assert(src_vec);
255         assert(tgt_vec);
256
257         src_len = src_vec->len;
258         tgt_len = tgt_vec->len;
259         assert(src_len > 0);
260         assert(tgt_len > 0);
261
262         mat = edge->costs;
263         assert(mat);
264
265         if (src_len == 1 && tgt_len == 1) {
266                 //panic("Something is wrong");
267         }
268
269         int allRowsOk = 1;
270         for (src_index = 0; src_index < src_len; ++src_index) {
271                 int onlyOneZero = 0;
272                 if (src_vec->entries[src_index].data == INF_COSTS) {
273                         continue;
274                 }
275                 for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
276                         if (tgt_vec->entries[tgt_index].data == INF_COSTS) {
277                                 continue;
278                         }
279                         if (mat->entries[src_index * tgt_len + tgt_index] == 0) {
280                                 if (onlyOneZero) {
281                                         onlyOneZero = 0;
282                                         break;
283                                 } else {
284                                         onlyOneZero = 1;
285                                         continue;
286                                 }
287                         }
288                         if (mat->entries[src_index * tgt_len + tgt_index] == INF_COSTS) {
289                                 continue;
290                         }
291                         onlyOneZero = 0;
292                         break;
293                 }
294                 allRowsOk &= onlyOneZero;
295         }
296
297         int allColsOk = 1;
298         for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
299                 int onlyOneZero = 0;
300                 if (tgt_vec->entries[tgt_index].data == INF_COSTS) {
301                         continue;
302                 }
303                 for (src_index = 0; src_index < src_len; ++src_index) {
304                         if (src_vec->entries[src_index].data == INF_COSTS) {
305                                 continue;
306                         }
307                         if (mat->entries[src_index * tgt_len + tgt_index] == 0) {
308                                 if (onlyOneZero) {
309                                         onlyOneZero = 0;
310                                         break;
311                                 } else {
312                                         onlyOneZero = 1;
313                                         continue;
314                                 }
315                         }
316                         if (mat->entries[src_index * tgt_len + tgt_index] == INF_COSTS) {
317                                 continue;
318                         }
319                         onlyOneZero = 0;
320                         break;
321                 }
322                 allColsOk &= onlyOneZero;
323         }
324
325         if (allRowsOk && allColsOk) {
326                 panic("Hurray");
327         }
328 }
329
330 static void simplify_edge(pbqp *pbqp, pbqp_edge *edge)
331 {
332         pbqp_matrix    *mat;
333         pbqp_node      *src_node;
334         pbqp_node      *tgt_node;
335         vector         *src_vec;
336         vector         *tgt_vec;
337         int             src_len;
338         int             tgt_len;
339
340         assert(pbqp);
341         assert(edge);
342
343         src_node = edge->src;
344         tgt_node = edge->tgt;
345         assert(src_node);
346         assert(tgt_node);
347
348         /* If edge are already deleted, we have nothing to do. */
349         if (!is_connected(src_node, edge) || !is_connected(tgt_node, edge))
350                 return;
351
352         if (pbqp->dump_file) {
353                 char txt[100];
354                 sprintf(txt, "Simplification of Edge n%d-n%d", src_node->index, tgt_node->index);
355                 dump_section(pbqp->dump_file, 3, txt);
356         }
357
358         src_vec = src_node->costs;
359         tgt_vec = tgt_node->costs;
360         assert(src_vec);
361         assert(tgt_vec);
362
363         src_len = src_vec->len;
364         tgt_len = tgt_vec->len;
365         assert(src_len > 0);
366         assert(tgt_len > 0);
367
368         mat = edge->costs;
369         assert(mat);
370
371         if (pbqp->dump_file) {
372                 fputs("Input:<br>\n", pbqp->dump_file);
373                 dump_simplifyedge(pbqp, edge);
374         }
375
376         normalize_towards_source(pbqp, edge);
377         normalize_towards_target(pbqp, edge);
378
379         if (pbqp->dump_file) {
380                 fputs("<br>\nOutput:<br>\n", pbqp->dump_file);
381                 dump_simplifyedge(pbqp, edge);
382         }
383
384         if (pbqp_matrix_is_zero(mat, src_vec, tgt_vec)) {
385                 if (pbqp->dump_file) {
386                         fputs("edge has been eliminated<br>\n", pbqp->dump_file);
387                 }
388
389                 delete_edge(edge);
390                 reorder_node(src_node);
391                 reorder_node(tgt_node);
392         } else {
393                 //check_melting_possibility(pbqp, edge);
394         }
395 }
396
397 static void initial_simplify_edges(pbqp *pbqp)
398 {
399         unsigned node_index;
400         unsigned node_len;
401
402         assert(pbqp);
403
404         if (pbqp->dump_file) {
405                 pbqp_dump_input(pbqp);
406                 dump_section(pbqp->dump_file, 1, "2. Simplification of Cost Matrices");
407         }
408
409         node_len = pbqp->num_nodes;
410
411         init_buckets();
412
413         /* First simplify all edges. */
414         for (node_index = 0; node_index < node_len; ++node_index) {
415                 unsigned    edge_index;
416                 pbqp_node  *node = get_node(pbqp, node_index);
417                 pbqp_edge **edges;
418                 unsigned    edge_len;
419
420                 if (!node) continue;
421
422                 edges = node->edges;
423                 edge_len = pbqp_node_get_degree(node);
424
425                 for (edge_index = 0; edge_index < edge_len; ++edge_index) {
426                         pbqp_edge *edge = edges[edge_index];
427
428                         /* Simplify only once per edge. */
429                         if (node != edge->src) continue;
430
431                         simplify_edge(pbqp, edge);
432                 }
433         }
434 }
435
436 num determine_solution(FILE *file)
437 {
438         unsigned node_index;
439         unsigned node_len;
440         num      solution   = 0;
441
442         if (file) {
443                 dump_section(file, 1, "4. Determine Solution/Minimum");
444                 dump_section(file, 2, "4.1. Trivial Solution");
445         }
446
447         /* Solve trivial nodes and calculate solution. */
448         node_len = node_bucket_get_length(node_buckets[0]);
449         for (node_index = 0; node_index < node_len; ++node_index) {
450                 pbqp_node *node = node_buckets[0][node_index];
451                 assert(node);
452
453                 node->solution = vector_get_min_index(node->costs);
454                 solution       = pbqp_add(solution,
455                                 node->costs->entries[node->solution].data);
456                 if (file) {
457                         fprintf(file, "node n%d is set to %d<br>\n", node->index, node->solution);
458                         dump_node(file, node);
459                 }
460         }
461
462         if (file) {
463                 dump_section(file, 2, "Minimum");
464                 fprintf(file, "Minimum is equal to %lld.", solution);
465         }
466
467         return solution;
468 }
469
470 static void back_propagate(pbqp *pbqp)
471 {
472         unsigned node_index;
473         unsigned node_len   = node_bucket_get_length(reduced_bucket);
474
475         assert(pbqp);
476         if (pbqp->dump_file) {
477                 dump_section(pbqp->dump_file, 2, "Back Propagation");
478         }
479
480         for (node_index = node_len; node_index > 0; --node_index) {
481                 pbqp_node *node = reduced_bucket[node_index - 1];
482
483                 switch (pbqp_node_get_degree(node)) {
484                         case 1:
485                                 back_propagate_RI(pbqp, node);
486                                 break;
487                         case 2:
488                                 back_propagate_RII(pbqp, node);
489                                 break;
490                         default:
491                                 panic("Only nodes with degree one or two should be in this bucket");
492                                 break;
493                 }
494         }
495 }
496
497 static void apply_heuristic_reductions(pbqp *pbqp)
498 {
499         for (;;) {
500                 if (edge_bucket_get_length(edge_bucket) > 0) {
501                         apply_edge(pbqp);
502                 } else if (node_bucket_get_length(node_buckets[1]) > 0) {
503                         apply_RI(pbqp);
504                 } else if (node_bucket_get_length(node_buckets[2]) > 0) {
505                         apply_RII(pbqp);
506                 } else if (node_bucket_get_length(node_buckets[3]) > 0) {
507                         apply_RN(pbqp);
508                 } else {
509                         return;
510                 }
511         }
512 }
513
514 void solve_pbqp_heuristical(pbqp *pbqp)
515 {
516         /* Reduce nodes degree ... */
517         initial_simplify_edges(pbqp);
518
519         /* ... and put node into bucket representing their degree. */
520         fill_node_buckets(pbqp);
521
522         FILE *fh = fopen("solutions.pb", "a");
523         fprintf(fh, "Solution");
524         fclose(fh);
525
526         apply_heuristic_reductions(pbqp);
527
528         pbqp->solution = determine_solution(pbqp->dump_file);
529
530         fh = fopen("solutions.pb", "a");
531         fprintf(fh, ": %lld\n", pbqp->solution);
532         fclose(fh);
533
534         /* Solve reduced nodes. */
535         back_propagate(pbqp);
536
537         free_buckets();
538 }
539
540 void apply_edge(pbqp *pbqp)
541 {
542         pbqp_edge *edge = edge_bucket_pop(&edge_bucket);
543
544         simplify_edge(pbqp, edge);
545 }
546
547 void apply_RI(pbqp *pbqp)
548 {
549         pbqp_node   *node       = node_bucket_pop(&node_buckets[1]);
550         pbqp_edge   *edge       = node->edges[0];
551         pbqp_matrix *mat        = edge->costs;
552         int          is_src     = edge->src == node;
553         pbqp_node   *other_node;
554
555         assert(pbqp_node_get_degree(node) == 1);
556
557         if (is_src) {
558                 other_node = edge->tgt;
559         } else {
560                 other_node = edge->src;
561         }
562
563         if (pbqp->dump_file) {
564                 char     txt[100];
565                 sprintf(txt, "RI-Reduction of Node n%d", node->index);
566                 dump_section(pbqp->dump_file, 2, txt);
567                 pbqp_dump_graph(pbqp);
568                 fputs("<br>\nBefore reduction:<br>\n", pbqp->dump_file);
569                 dump_node(pbqp->dump_file, node);
570                 dump_node(pbqp->dump_file, other_node);
571                 dump_edge(pbqp->dump_file, edge);
572         }
573
574         if (is_src) {
575                 pbqp_matrix_add_to_all_cols(mat, node->costs);
576                 normalize_towards_target(pbqp, edge);
577         } else {
578                 pbqp_matrix_add_to_all_rows(mat, node->costs);
579                 normalize_towards_source(pbqp, edge);
580         }
581         disconnect_edge(other_node, edge);
582
583         if (pbqp->dump_file) {
584                 fputs("<br>\nAfter reduction:<br>\n", pbqp->dump_file);
585                 dump_node(pbqp->dump_file, other_node);
586         }
587
588         reorder_node(other_node);
589
590         /* Add node to back propagation list. */
591         node_bucket_insert(&reduced_bucket, node);
592 }
593
594 void apply_RII(pbqp *pbqp)
595 {
596         pbqp_node   *node       = node_bucket_pop(&node_buckets[2]);
597         pbqp_edge   *src_edge   = node->edges[0];
598         pbqp_edge   *tgt_edge   = node->edges[1];
599         int          src_is_src = src_edge->src == node;
600         int          tgt_is_src = tgt_edge->src == node;
601         pbqp_matrix *src_mat;
602         pbqp_matrix *tgt_mat;
603         pbqp_node   *src_node;
604         pbqp_node   *tgt_node;
605         pbqp_matrix *mat;
606         vector      *vec;
607         vector      *node_vec;
608         vector      *src_vec;
609         vector      *tgt_vec;
610         unsigned     col_index;
611         unsigned     col_len;
612         unsigned     row_index;
613         unsigned     row_len;
614         unsigned     node_len;
615
616         assert(pbqp);
617         assert(pbqp_node_get_degree(node) == 2);
618
619         if (src_is_src) {
620                 src_node = src_edge->tgt;
621         } else {
622                 src_node = src_edge->src;
623         }
624
625         if (tgt_is_src) {
626                 tgt_node = tgt_edge->tgt;
627         } else {
628                 tgt_node = tgt_edge->src;
629         }
630
631         /* Swap nodes if necessary. */
632         if (tgt_node->index < src_node->index) {
633                 pbqp_node *tmp_node;
634                 pbqp_edge *tmp_edge;
635
636                 tmp_node = src_node;
637                 src_node = tgt_node;
638                 tgt_node = tmp_node;
639
640                 tmp_edge = src_edge;
641                 src_edge = tgt_edge;
642                 tgt_edge = tmp_edge;
643
644                 src_is_src = src_edge->src == node;
645                 tgt_is_src = tgt_edge->src == node;
646         }
647
648         if (pbqp->dump_file) {
649                 char     txt[100];
650                 sprintf(txt, "RII-Reduction of Node n%d", node->index);
651                 dump_section(pbqp->dump_file, 2, txt);
652                 pbqp_dump_graph(pbqp);
653                 fputs("<br>\nBefore reduction:<br>\n", pbqp->dump_file);
654                 dump_node(pbqp->dump_file, src_node);
655                 dump_edge(pbqp->dump_file, src_edge);
656                 dump_node(pbqp->dump_file, node);
657                 dump_edge(pbqp->dump_file, tgt_edge);
658                 dump_node(pbqp->dump_file, tgt_node);
659         }
660
661         src_mat = src_edge->costs;
662         tgt_mat = tgt_edge->costs;
663
664         src_vec  = src_node->costs;
665         tgt_vec  = tgt_node->costs;
666         node_vec = node->costs;
667
668         row_len  = src_vec->len;
669         col_len  = tgt_vec->len;
670         node_len = node_vec->len;
671
672         mat = pbqp_matrix_alloc(pbqp, row_len, col_len);
673
674         for (row_index = 0; row_index < row_len; ++row_index) {
675                 for (col_index = 0; col_index < col_len; ++col_index) {
676                         vec = vector_copy(pbqp, node_vec);
677
678                         if (src_is_src) {
679                                 vector_add_matrix_col(vec, src_mat, row_index);
680                         } else {
681                                 vector_add_matrix_row(vec, src_mat, row_index);
682                         }
683
684                         if (tgt_is_src) {
685                                 vector_add_matrix_col(vec, tgt_mat, col_index);
686                         } else {
687                                 vector_add_matrix_row(vec, tgt_mat, col_index);
688                         }
689
690                         mat->entries[row_index * col_len + col_index] = vector_get_min(vec);
691
692                         obstack_free(&pbqp->obstack, vec);
693                 }
694         }
695
696         pbqp_edge *edge = get_edge(pbqp, src_node->index, tgt_node->index);
697
698         /* Disconnect node. */
699         disconnect_edge(src_node, src_edge);
700         disconnect_edge(tgt_node, tgt_edge);
701
702         /* Add node to back propagation list. */
703         node_bucket_insert(&reduced_bucket, node);
704
705         if (edge == NULL) {
706                 edge = alloc_edge(pbqp, src_node->index, tgt_node->index, mat);
707         } else {
708                 pbqp_matrix_add(edge->costs, mat);
709
710                 /* Free local matrix. */
711                 obstack_free(&pbqp->obstack, mat);
712
713                 reorder_node(src_node);
714                 reorder_node(tgt_node);
715         }
716
717         if (pbqp->dump_file) {
718                 fputs("<br>\nAfter reduction:<br>\n", pbqp->dump_file);
719                 dump_edge(pbqp->dump_file, edge);
720         }
721
722         /* Edge has changed so we simplify it. */
723         simplify_edge(pbqp, edge);
724 }
725
726 static void select_alternative(pbqp_node *node, unsigned selected_index)
727 {
728         unsigned  edge_index;
729         unsigned  node_index;
730         unsigned  node_len;
731         vector   *node_vec;
732         unsigned  max_degree = pbqp_node_get_degree(node);
733
734         assert(node);
735         node->solution = selected_index;
736         node_vec = node->costs;
737         node_len = node_vec->len;
738         assert(selected_index < node_len);
739
740         /* Set all other costs to infinity. */
741         for (node_index = 0; node_index < node_len; ++node_index) {
742                 if (node_index != selected_index) {
743                         node_vec->entries[node_index].data = INF_COSTS;
744                 }
745         }
746
747         /* Add all incident edges to edge bucket, since they are now independent. */
748         for (edge_index = 0; edge_index < max_degree; ++edge_index) {
749                 insert_into_edge_bucket(node->edges[edge_index]);
750         }
751 }
752
753 static pbqp_node *get_node_with_max_degree(void)
754 {
755         pbqp_node  **bucket       = node_buckets[3];
756         unsigned     bucket_len   = node_bucket_get_length(bucket);
757         unsigned     bucket_index;
758         unsigned     max_degree   = 0;
759         pbqp_node   *result       = NULL;
760
761         for (bucket_index = 0; bucket_index < bucket_len; ++bucket_index) {
762                 pbqp_node *candidate = bucket[bucket_index];
763                 unsigned   degree    = pbqp_node_get_degree(candidate);
764
765                 if (degree > max_degree) {
766                         result = candidate;
767                         max_degree = degree;
768                 }
769         }
770
771         return result;
772 }
773
774 static unsigned get_local_minimal_alternative(pbqp *pbqp, pbqp_node *node)
775 {
776         pbqp_edge   *edge;
777         vector      *node_vec;
778         vector      *vec;
779         pbqp_matrix *mat;
780         unsigned     edge_index;
781         unsigned     max_degree   = 0;
782         unsigned     node_index;
783         unsigned     node_len;
784         unsigned     min_index    = 0;
785         num          min          = INF_COSTS;
786         int          is_src;
787
788         assert(pbqp);
789         assert(node);
790         node_vec = node->costs;
791         node_len = node_vec->len;
792
793         for (node_index = 0; node_index < node_len; ++node_index) {
794                 num value = node_vec->entries[node_index].data;
795
796                 for (edge_index = 0; edge_index < max_degree; ++edge_index) {
797                         edge   = node->edges[edge_index];
798                         mat    = edge->costs;
799                         is_src = edge->src == node;
800
801                         if (is_src) {
802                                 vec = vector_copy(pbqp, edge->tgt->costs);
803                                 vector_add_matrix_row(vec, mat, node_index);
804                         } else {
805                                 vec = vector_copy(pbqp, edge->src->costs);
806                                 vector_add_matrix_col(vec, mat, node_index);
807                         }
808
809                         value = pbqp_add(value, vector_get_min(vec));
810
811                         obstack_free(&pbqp->obstack, vec);
812                 }
813
814                 if (value < min) {
815                         min = value;
816                         min_index = node_index;
817                 }
818         }
819
820         return min_index;
821 }
822
823 void apply_RN(pbqp *pbqp)
824 {
825         pbqp_node   *node         = NULL;
826         unsigned     min_index    = 0;
827
828         assert(pbqp);
829
830         /* We want to reduce a node with maximum degree. */
831         node = get_node_with_max_degree();
832         assert(node);
833         assert(pbqp_node_get_degree(node) > 2);
834
835         if (pbqp->dump_file) {
836                 char     txt[100];
837                 sprintf(txt, "RN-Reduction of Node n%d", node->index);
838                 dump_section(pbqp->dump_file, 2, txt);
839                 pbqp_dump_graph(pbqp);
840         }
841
842         min_index = get_local_minimal_alternative(pbqp, node);
843
844         if (pbqp->dump_file) {
845                 fprintf(pbqp->dump_file, "node n%d is set to %d<br><br>\n",
846                                         node->index, min_index);
847         }
848
849         FILE *fh = fopen("solutions.pb", "a");
850         fprintf(fh, "[%u]", min_index);
851         fclose(fh);
852
853         /* Now that we found the local minimum set all other costs to infinity. */
854         select_alternative(node, min_index);
855 }
856
857 static void apply_brute_force_reductions(pbqp *pbqp)
858 {
859         for (;;) {
860                 if (edge_bucket_get_length(edge_bucket) > 0) {
861                         apply_edge(pbqp);
862                 } else if (node_bucket_get_length(node_buckets[1]) > 0) {
863                         apply_RI(pbqp);
864                 } else if (node_bucket_get_length(node_buckets[2]) > 0) {
865                         apply_RII(pbqp);
866                 } else if (node_bucket_get_length(node_buckets[3]) > 0) {
867                         apply_Brute_Force(pbqp);
868                 } else {
869                         return;
870                 }
871         }
872 }
873
874 static unsigned get_minimal_alternative(pbqp *pbqp, pbqp_node *node)
875 {
876         vector      *node_vec;
877         unsigned     node_index;
878         unsigned     node_len;
879         unsigned     min_index    = 0;
880         num          min          = INF_COSTS;
881         unsigned     bucket_index;
882
883         assert(pbqp);
884         assert(node);
885         node_vec     = node->costs;
886         node_len     = node_vec->len;
887         bucket_index = node->bucket_index;
888
889         for (node_index = 0; node_index < node_len; ++node_index) {
890                 pbqp_node_bucket bucket_deg3;
891                 num              value;
892                 unsigned         bucket_0_length;
893                 unsigned         bucket_red_length;
894
895                 char *tmp = obstack_finish(&pbqp->obstack);
896
897                 node_bucket_init(&bucket_deg3);
898
899                 /* Some node buckets and the edge bucket should be empty. */
900                 assert(node_bucket_get_length(node_buckets[1]) == 0);
901                 assert(node_bucket_get_length(node_buckets[2]) == 0);
902                 assert(edge_bucket_get_length(edge_bucket)     == 0);
903
904                 /* char *tmp = obstack_finish(&pbqp->obstack); */
905
906                 /* Save current PBQP state. */
907                 node_bucket_copy(&bucket_deg3, node_buckets[3]);
908                 node_bucket_shrink(&node_buckets[3], 0);
909                 node_bucket_deep_copy(pbqp, &node_buckets[3], bucket_deg3);
910                 node_bucket_update(pbqp, node_buckets[3]);
911                 bucket_0_length   = node_bucket_get_length(node_buckets[0]);
912                 bucket_red_length = node_bucket_get_length(reduced_bucket);
913
914                 /* Select alternative and solve PBQP recursively. */
915                 select_alternative(node_buckets[3][bucket_index], node_index);
916                 apply_brute_force_reductions(pbqp);
917
918                 value = determine_solution(pbqp->dump_file);
919
920                 if (value < min) {
921                         min = value;
922                         min_index = node_index;
923                 }
924
925                 /* Some node buckets and the edge bucket should still be empty. */
926                 assert(node_bucket_get_length(node_buckets[1]) == 0);
927                 assert(node_bucket_get_length(node_buckets[2]) == 0);
928                 assert(edge_bucket_get_length(edge_bucket)     == 0);
929
930                 /* Clear modified buckets... */
931                 node_bucket_shrink(&node_buckets[3], 0);
932
933                 /* ... and restore old PBQP state. */
934                 node_bucket_shrink(&node_buckets[0], bucket_0_length);
935                 node_bucket_shrink(&reduced_bucket, bucket_red_length);
936                 node_bucket_copy(&node_buckets[3], bucket_deg3);
937                 node_bucket_update(pbqp, node_buckets[3]);
938
939                 /* Free copies. */
940                 /* obstack_free(&pbqp->obstack, tmp); */
941                 node_bucket_free(&bucket_deg3);
942
943                 obstack_free(&pbqp->obstack, tmp);
944         }
945
946         return min_index;
947 }
948
949 void apply_Brute_Force(pbqp *pbqp)
950 {
951         pbqp_node   *node         = NULL;
952         unsigned     min_index    = 0;
953
954         assert(pbqp);
955
956         /* We want to reduce a node with maximum degree. */
957         node = get_node_with_max_degree();
958         assert(node);
959         assert(pbqp_node_get_degree(node) > 2);
960
961         if (pbqp->dump_file) {
962                 char     txt[100];
963                 sprintf(txt, "BF-Reduction of Node n%d", node->index);
964                 dump_section(pbqp->dump_file, 2, txt);
965                 pbqp_dump_graph(pbqp);
966         }
967
968         dump++;
969         min_index = get_minimal_alternative(pbqp, node);
970         node = pbqp->nodes[node->index];
971
972         if (pbqp->dump_file) {
973                 fprintf(pbqp->dump_file, "node n%d is set to %d<br><br>\n",
974                                         node->index, min_index);
975         }
976
977         dump--;
978         if (dump == 0) {
979                 FILE *fh = fopen("solutions.pb", "a");
980                 fprintf(fh, "[%u]", min_index);
981                 fclose(fh);
982         }
983
984         /* Now that we found the minimum set all other costs to infinity. */
985         select_alternative(node, min_index);
986 }
987
988 void solve_pbqp_brute_force(pbqp *pbqp)
989 {
990         /* Reduce nodes degree ... */
991         initial_simplify_edges(pbqp);
992
993         /* ... and put node into bucket representing their degree. */
994         fill_node_buckets(pbqp);
995
996         FILE *fh = fopen("solutions.pb", "a");
997         fprintf(fh, "Solution");
998         fclose(fh);
999
1000         apply_brute_force_reductions(pbqp);
1001
1002         pbqp->solution = determine_solution(pbqp->dump_file);
1003
1004         fh = fopen("solutions.pb", "a");
1005         fprintf(fh, ": %lld\n", pbqp->solution);
1006         fclose(fh);
1007
1008         /* Solve reduced nodes. */
1009         back_propagate(pbqp);
1010
1011         free_buckets();
1012 }
1013
1014 void back_propagate_RI(pbqp *pbqp, pbqp_node *node)
1015 {
1016         pbqp_edge   *edge;
1017         pbqp_node   *other;
1018         pbqp_matrix *mat;
1019         vector      *vec;
1020         int          is_src;
1021
1022         assert(pbqp);
1023         assert(node);
1024
1025         edge = node->edges[0];
1026         mat = edge->costs;
1027         is_src = edge->src == node;
1028         vec = node->costs;
1029
1030         if (is_src) {
1031                 other = edge->tgt;
1032                 assert(other);
1033
1034                 /* Update pointer for brute force solver. */
1035                 other = pbqp->nodes[other->index];
1036
1037                 vector_add_matrix_col(vec, mat, other->solution);
1038         } else {
1039                 other = edge->src;
1040                 assert(other);
1041
1042                 /* Update pointer for brute force solver. */
1043                 other = pbqp->nodes[other->index];
1044
1045                 vector_add_matrix_row(vec, mat, other->solution);
1046         }
1047
1048         node->solution = vector_get_min_index(vec);
1049         if (pbqp->dump_file) {
1050                 fprintf(pbqp->dump_file, "node n%d is set to %d<br>\n", node->index, node->solution);
1051         }
1052 }
1053
1054 void back_propagate_RII(pbqp *pbqp, pbqp_node *node)
1055 {
1056         pbqp_edge   *src_edge   = node->edges[0];
1057         pbqp_edge   *tgt_edge   = node->edges[1];
1058         int          src_is_src = src_edge->src == node;
1059         int          tgt_is_src = tgt_edge->src == node;
1060         pbqp_matrix *src_mat;
1061         pbqp_matrix *tgt_mat;
1062         pbqp_node   *src_node;
1063         pbqp_node   *tgt_node;
1064         vector      *vec;
1065         vector      *node_vec;
1066         unsigned     col_index;
1067         unsigned     row_index;
1068
1069         assert(pbqp);
1070
1071         if (src_is_src) {
1072                 src_node = src_edge->tgt;
1073         } else {
1074                 src_node = src_edge->src;
1075         }
1076
1077         if (tgt_is_src) {
1078                 tgt_node = tgt_edge->tgt;
1079         } else {
1080                 tgt_node = tgt_edge->src;
1081         }
1082
1083         /* Swap nodes if necessary. */
1084         if (tgt_node->index < src_node->index) {
1085                 pbqp_node *tmp_node;
1086                 pbqp_edge *tmp_edge;
1087
1088                 tmp_node = src_node;
1089                 src_node = tgt_node;
1090                 tgt_node = tmp_node;
1091
1092                 tmp_edge = src_edge;
1093                 src_edge = tgt_edge;
1094                 tgt_edge = tmp_edge;
1095
1096                 src_is_src = src_edge->src == node;
1097                 tgt_is_src = tgt_edge->src == node;
1098         }
1099
1100         /* Update pointer for brute force solver. */
1101         src_node = pbqp->nodes[src_node->index];
1102         tgt_node = pbqp->nodes[tgt_node->index];
1103
1104         src_mat = src_edge->costs;
1105         tgt_mat = tgt_edge->costs;
1106
1107         node_vec = node->costs;
1108
1109         row_index = src_node->solution;
1110         col_index = tgt_node->solution;
1111
1112         vec = vector_copy(pbqp, node_vec);
1113
1114         if (src_is_src) {
1115                 vector_add_matrix_col(vec, src_mat, row_index);
1116         } else {
1117                 vector_add_matrix_row(vec, src_mat, row_index);
1118         }
1119
1120         if (tgt_is_src) {
1121                 vector_add_matrix_col(vec, tgt_mat, col_index);
1122         } else {
1123                 vector_add_matrix_row(vec, tgt_mat, col_index);
1124         }
1125
1126         node->solution = vector_get_min_index(vec);
1127         if (pbqp->dump_file) {
1128                 fprintf(pbqp->dump_file, "node n%d is set to %d<br>\n", node->index, node->solution);
1129         }
1130
1131         obstack_free(&pbqp->obstack, vec);
1132 }
1133
1134 int node_is_reduced(pbqp_node *node)
1135 {
1136         if (!reduced_bucket) return 0;
1137
1138         if (pbqp_node_get_degree(node) == 0) return 1;
1139
1140         return node_bucket_contains(reduced_bucket, node);
1141 }