5 #include "heuristical.h"
6 #include "html_dumper.h"
10 #include "pbqp_edge_t.h"
11 #include "pbqp_node.h"
12 #include "pbqp_node_t.h"
15 static pbqp_edge **edge_bucket;
16 static pbqp_node **node_buckets[4];
17 static pbqp_node **reduced_bucket = NULL;
18 static int buckets_filled = 0;
20 static void init_buckets(void)
24 edge_bucket = NEW_ARR_F(pbqp_edge *, 0);
25 reduced_bucket = NEW_ARR_F(pbqp_node *, 0);
27 for (i = 0; i < 4; ++i) {
28 node_buckets[i] = NEW_ARR_F(pbqp_node *, 0);
32 static void free_buckets(void)
36 for (i = 0; i < 4; ++i) {
37 DEL_ARR_F(node_buckets[i]);
38 node_buckets[i] = NULL;
41 DEL_ARR_F(edge_bucket);
44 DEL_ARR_F(reduced_bucket);
45 reduced_bucket = NULL;
50 static void fill_node_buckets(pbqp *pbqp)
56 node_len = pbqp->num_nodes;
58 for (node_index = 0; node_index < node_len; ++node_index) {
60 pbqp_node *node = get_node(pbqp, node_index);
64 arity = ARR_LEN(node->edges);
66 /* We have only one bucket for nodes with arity >= 3. */
71 node->bucket_index = ARR_LEN(node_buckets[arity]);
73 ARR_APP1(pbqp_node *, node_buckets[arity], node);
79 static void normalize_towards_source(pbqp *pbqp, pbqp_edge *edge)
98 src_vec = src_node->costs;
99 tgt_vec = tgt_node->costs;
103 src_len = src_vec->len;
104 tgt_len = tgt_vec->len;
111 /* Normalize towards source node. */
112 for (src_index = 0; src_index < src_len; ++src_index) {
113 num min = pbqp_matrix_get_row_min(mat, src_index, tgt_vec);
116 pbqp_matrix_sub_row_value(mat, src_index, tgt_vec, min);
117 src_vec->entries[src_index].data = pbqp_add(
118 src_vec->entries[src_index].data, min);
120 // TODO add to edge_list if inf
125 static void normalize_towards_target(pbqp *pbqp, pbqp_edge *edge)
139 src_node = edge->src;
140 tgt_node = edge->tgt;
144 src_vec = src_node->costs;
145 tgt_vec = tgt_node->costs;
149 src_len = src_vec->len;
150 tgt_len = tgt_vec->len;
157 for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
158 num min = pbqp_matrix_get_col_min(mat, tgt_index, src_vec);
161 pbqp_matrix_sub_col_value(mat, tgt_index, src_vec, min);
162 tgt_vec->entries[tgt_index].data = pbqp_add(
163 tgt_vec->entries[tgt_index].data, min);
165 // TODO add to edge_list if inf
170 static void reorder_node(pbqp_node *node)
174 unsigned old_bucket_len;
175 unsigned old_bucket_index;
176 pbqp_node **old_bucket;
179 if (!buckets_filled) return;
183 arity = ARR_LEN(node->edges);
185 /* Same bucket as before */
186 if (arity > 2) return;
188 /* Assume node lost one incident edge. */
189 old_arity = arity + 1;
190 old_bucket = node_buckets[old_arity];
191 old_bucket_len = ARR_LEN(old_bucket);
192 old_bucket_index = node->bucket_index;
194 if (old_bucket_len <= old_bucket_index ||
195 old_bucket[old_bucket_index] != node) {
196 /* Old arity is new arity, so we have nothing to do. */
197 assert(old_bucket_index < ARR_LEN(node_buckets[arity]) &&
198 node_buckets[arity][old_bucket_index] == node);
202 assert(old_bucket[old_bucket_index] == node);
204 /* Delete node from old bucket... */
205 other = old_bucket[old_bucket_len - 1];
206 other->bucket_index = old_bucket_index;
207 old_bucket[old_bucket_index] = other;
208 ARR_SHRINKLEN(node_buckets[old_arity], old_bucket_len - 1);
210 /* ..and add to new one. */
211 node->bucket_index = ARR_LEN(node_buckets[arity]);
212 ARR_APP1(pbqp_node*, node_buckets[arity], node);
215 static void simplify_edge(pbqp *pbqp, pbqp_edge *edge)
228 src_node = edge->src;
229 tgt_node = edge->tgt;
233 if (pbqp->dump_file) {
235 sprintf(txt, "Simplification of Edge n%d-n%d", src_node->index, tgt_node->index);
236 dump_section(pbqp->dump_file, 3, txt);
239 src_vec = src_node->costs;
240 tgt_vec = tgt_node->costs;
244 src_len = src_vec->len;
245 tgt_len = tgt_vec->len;
252 if (pbqp->dump_file) {
253 fputs("Input:<br>\n", pbqp->dump_file);
254 dump_simplifyedge(pbqp, edge);
257 normalize_towards_source(pbqp, edge);
258 normalize_towards_target(pbqp, edge);
260 if (pbqp->dump_file) {
261 fputs("<br>\nOutput:<br>\n", pbqp->dump_file);
262 dump_simplifyedge(pbqp, edge);
265 if (pbqp_matrix_is_zero(mat, src_vec, tgt_vec)) {
266 if (pbqp->dump_file) {
267 fputs("edge has been eliminated<br>\n", pbqp->dump_file);
271 reorder_node(src_node);
272 reorder_node(tgt_node);
276 void solve_pbqp_heuristical(pbqp *pbqp)
283 if (pbqp->dump_file) {
284 pbqp_dump_input(pbqp);
285 dump_section(pbqp->dump_file, 1, "2. Simplification of Cost Matrices");
288 node_len = pbqp->num_nodes;
292 /* First simplify all edges. */
293 for (node_index = 0; node_index < node_len; ++node_index) {
295 pbqp_node *node = get_node(pbqp, node_index);
302 edge_len = ARR_LEN(edges);
304 for (edge_index = 0; edge_index < edge_len; ++edge_index) {
305 pbqp_edge *edge = edges[edge_index];
307 /* Simplify only once per edge. */
308 if (node != edge->src) continue;
310 simplify_edge(pbqp, edge);
314 /* Put node into bucket representing their arity. */
315 fill_node_buckets(pbqp);
318 if (ARR_LEN(edge_bucket) > 0) {
319 panic("Please implement edge simplification");
320 } else if (ARR_LEN(node_buckets[1]) > 0) {
322 } else if (ARR_LEN(node_buckets[2]) > 0) {
324 } else if (ARR_LEN(node_buckets[3]) > 0) {
325 panic("Please implement RN simplification");
331 if (pbqp->dump_file) {
332 dump_section(pbqp->dump_file, 1, "4. Determine Solution/Minimum");
333 dump_section(pbqp->dump_file, 2, "4.1. Trivial Solution");
336 /* Solve trivial nodes and calculate solution. */
337 node_len = ARR_LEN(node_buckets[0]);
338 for (node_index = 0; node_index < node_len; ++node_index) {
339 pbqp_node *node = node_buckets[0][node_index];
342 node->solution = vector_get_min_index(node->costs);
343 pbqp->solution = pbqp_add(pbqp->solution,
344 node->costs->entries[node->solution].data);
345 if (pbqp->dump_file) {
346 fprintf(pbqp->dump_file, "node n%d is set to %d<br>\n", node->index, node->solution);
347 dump_node(pbqp, node);
351 if (pbqp->dump_file) {
352 dump_section(pbqp->dump_file, 2, "Minimum");
353 fprintf(pbqp->dump_file, "Minimum is equal to %d.", pbqp->solution);
354 dump_section(pbqp->dump_file, 2, "Back Propagation");
357 /* Solve reduced nodes. */
358 node_len = ARR_LEN(reduced_bucket);
359 for (node_index = node_len; node_index > 0; --node_index) {
360 pbqp_node *node = reduced_bucket[node_index - 1];
363 switch (ARR_LEN(node->edges)) {
365 back_propagate_RI(pbqp, node);
368 back_propagate_RII(pbqp, node);
371 panic("Only nodes with degree one or two should be in this bucket");
379 void apply_RI(pbqp *pbqp)
381 pbqp_node **bucket = node_buckets[1];
382 unsigned bucket_len = ARR_LEN(bucket);
383 pbqp_node *node = bucket[bucket_len - 1];
384 pbqp_edge *edge = node->edges[0];
385 pbqp_matrix *mat = edge->costs;
386 int is_src = edge->src == node;
387 pbqp_node *other_node;
390 other_node = edge->tgt;
392 other_node = edge->src;
395 if (pbqp->dump_file) {
397 sprintf(txt, "RI-Reduktion of Node n%d", node->index);
398 dump_section(pbqp->dump_file, 2, txt);
399 pbqp_dump_graph(pbqp);
400 fputs("<br>\nBefore reduction:<br>\n", pbqp->dump_file);
401 dump_node(pbqp, node);
402 dump_node(pbqp, other_node);
403 dump_edge(pbqp, edge);
407 pbqp_matrix_add_to_all_cols(mat, node->costs);
408 normalize_towards_target(pbqp, edge);
410 pbqp_matrix_add_to_all_rows(mat, node->costs);
411 normalize_towards_source(pbqp, edge);
413 disconnect_edge(other_node, edge);
415 if (pbqp->dump_file) {
416 fputs("<br>\nAfter reduction:<br>\n", pbqp->dump_file);
417 dump_node(pbqp, other_node);
420 /* Remove node from bucket... */
421 ARR_SHRINKLEN(bucket, (int)bucket_len - 1);
422 reorder_node(other_node);
424 /* ...and add it to back propagation list. */
425 node->bucket_index = ARR_LEN(reduced_bucket);
426 ARR_APP1(pbqp_node *, reduced_bucket, node);
429 void apply_RII(pbqp *pbqp)
431 pbqp_node **bucket = node_buckets[2];
432 unsigned bucket_len = ARR_LEN(bucket);
433 pbqp_node *node = bucket[bucket_len - 1];
434 pbqp_edge *src_edge = node->edges[0];
435 pbqp_edge *tgt_edge = node->edges[1];
436 int src_is_src = src_edge->src == node;
437 int tgt_is_src = tgt_edge->src == node;
438 pbqp_matrix *src_mat;
439 pbqp_matrix *tgt_mat;
456 src_node = src_edge->tgt;
458 src_node = src_edge->src;
462 tgt_node = tgt_edge->tgt;
464 tgt_node = tgt_edge->src;
467 /* Swap nodes if necessary. */
468 if (tgt_node->index < src_node->index) {
480 src_is_src = src_edge->src == node;
481 tgt_is_src = tgt_edge->src == node;
484 if (pbqp->dump_file) {
486 sprintf(txt, "RII-Reduktion of Node n%d", node->index);
487 dump_section(pbqp->dump_file, 2, txt);
488 pbqp_dump_graph(pbqp);
489 fputs("<br>\nBefore reduction:<br>\n", pbqp->dump_file);
490 dump_node(pbqp, src_node);
491 dump_edge(pbqp, src_edge);
492 dump_node(pbqp, node);
493 dump_edge(pbqp, tgt_edge);
494 dump_node(pbqp, tgt_node);
497 src_mat = src_edge->costs;
498 tgt_mat = tgt_edge->costs;
500 src_vec = src_node->costs;
501 tgt_vec = tgt_node->costs;
502 node_vec = node->costs;
504 row_len = src_vec->len;
505 col_len = tgt_vec->len;
506 node_len = node_vec->len;
508 mat = pbqp_matrix_alloc(pbqp, row_len, col_len);
510 for (row_index = 0; row_index < row_len; ++row_index) {
511 for (col_index = 0; col_index < col_len; ++col_index) {
512 vec = vector_copy(pbqp, node_vec);
515 vector_add_matrix_col(vec, src_mat, row_index);
517 vector_add_matrix_row(vec, src_mat, row_index);
521 vector_add_matrix_col(vec, tgt_mat, col_index);
523 vector_add_matrix_row(vec, tgt_mat, col_index);
526 mat->entries[row_index * col_len + col_index] = vector_get_min(vec);
528 obstack_free(&pbqp->obstack, vec);
532 pbqp_edge *edge = get_edge(pbqp, src_node->index, tgt_node->index);
534 /* Disconnect node. */
535 disconnect_edge(src_node, src_edge);
536 disconnect_edge(tgt_node, tgt_edge);
538 /* Remove node from bucket... */
539 ARR_SHRINKLEN(bucket, (int)bucket_len - 1);
541 /* ...and add it to back propagation list. */
542 node->bucket_index = ARR_LEN(reduced_bucket);
543 ARR_APP1(pbqp_node *, reduced_bucket, node);
546 edge = alloc_edge(pbqp, src_node->index, tgt_node->index, mat);
548 pbqp_matrix_add(edge->costs, mat);
550 /* Free local matrix. */
551 obstack_free(&pbqp->obstack, mat);
553 reorder_node(src_node);
554 reorder_node(tgt_node);
557 if (pbqp->dump_file) {
558 fputs("<br>\nAfter reduction:<br>\n", pbqp->dump_file);
559 dump_edge(pbqp, edge);
562 /* Edge has changed so we simplify it. */
563 simplify_edge(pbqp, edge);
566 void apply_RN(pbqp *pbqp)
568 pbqp_node **bucket = node_buckets[3];
569 unsigned bucket_len = ARR_LEN(bucket);
570 pbqp_node *node = bucket[bucket_len - 1];
572 vector *node_vec = node->costs;
576 unsigned edge_len = ARR_LEN(node->edges);
578 unsigned node_len = ARR_LEN(node_vec);
579 unsigned min_index = 0;
585 for (node_index = 0; node_index < node_len; ++node_index) {
588 for (edge_index = 0; edge_index < edge_len; ++edge_index) {
589 edge = node->edges[edge_index];
591 is_src = edge->src == node;
594 vec = vector_copy(pbqp, edge->tgt->costs);
595 vector_add_matrix_row(vec, mat, node_index);
597 vec = vector_copy(pbqp, edge->src->costs);
598 vector_add_matrix_col(vec, mat, node_index);
601 value = pbqp_add(value, vector_get_min(vec));
603 obstack_free(&pbqp->obstack, vec);
608 min_index = node_index;
612 node->solution = min_index;
614 /* Now that we found the local minimum set all other costs to infinity. */
615 for (node_index = 0; node_index < node_len; ++node_index) {
616 if (node_index != min_index) {
617 node_vec->entries[node_index].data = INF_COSTS;
621 /* Add all incident edges to edge bucket, since they are now independent. */
622 for (edge_index = 0; edge_index < edge_len; ++edge_index) {
623 ARR_APP1(pbqp_edge *, edge_bucket, node->edges[node_index]);
627 void back_propagate_RI(pbqp *pbqp, pbqp_node *node)
638 edge = node->edges[0];
640 is_src = edge->src == node;
646 vector_add_matrix_col(vec, mat, other->solution);
650 vector_add_matrix_row(vec, mat, other->solution);
653 node->solution = vector_get_min_index(vec);
654 if (pbqp->dump_file) {
655 fprintf(pbqp->dump_file, "node n%d is set to %d<br>\n", node->index, node->solution);
659 void back_propagate_RII(pbqp *pbqp, pbqp_node *node)
661 pbqp_edge *src_edge = node->edges[0];
662 pbqp_edge *tgt_edge = node->edges[1];
663 int src_is_src = src_edge->src == node;
664 int tgt_is_src = tgt_edge->src == node;
665 pbqp_matrix *src_mat;
666 pbqp_matrix *tgt_mat;
677 src_node = src_edge->tgt;
679 src_node = src_edge->src;
683 tgt_node = tgt_edge->tgt;
685 tgt_node = tgt_edge->src;
688 /* Swap nodes if necessary. */
689 if (tgt_node->index < src_node->index) {
701 src_is_src = src_edge->src == node;
702 tgt_is_src = tgt_edge->src == node;
705 src_mat = src_edge->costs;
706 tgt_mat = tgt_edge->costs;
708 node_vec = node->costs;
710 row_index = src_node->solution;
711 col_index = tgt_node->solution;
713 vec = vector_copy(pbqp, node_vec);
716 vector_add_matrix_col(vec, src_mat, row_index);
718 vector_add_matrix_row(vec, src_mat, row_index);
722 vector_add_matrix_col(vec, tgt_mat, col_index);
724 vector_add_matrix_row(vec, tgt_mat, col_index);
727 node->solution = vector_get_min_index(vec);
728 if (pbqp->dump_file) {
729 fprintf(pbqp->dump_file, "node n%d is set to %d<br>\n", node->index, node->solution);
732 obstack_free(&pbqp->obstack, vec);
735 int node_is_reduced(pbqp_node *node)
737 if (!reduced_bucket) return 0;
740 if (ARR_LEN(node->edges) == 0) return 1;
742 unsigned bucket_length = ARR_LEN(reduced_bucket);
743 unsigned bucket_index = node->bucket_index;
745 return bucket_index < bucket_length && reduced_bucket[bucket_index] == node;