+#include "kaps.h"
+#include "matrix.h"
+#include "pbqp_edge.h"
+#include "pbqp_edge_t.h"
+#include "pbqp_node.h"
+#include "pbqp_node_t.h"
+#include "vector.h"
+
+static pbqp_edge **edge_bucket;
+static pbqp_node **node_buckets[4];
+static pbqp_node **reduced_bucket = NULL;
+static int buckets_filled = 0;
+
+static num add(num x, num y)
+{
+ if (x == INF_COSTS || y == INF_COSTS) return INF_COSTS;
+
+ return x + y;
+}
+
+static void init_buckets(void)
+{
+ int i;
+
+ edge_bucket = NEW_ARR_F(pbqp_edge *, 0);
+ reduced_bucket = NEW_ARR_F(pbqp_node *, 0);
+
+ for (i = 0; i < 4; ++i) {
+ node_buckets[i] = NEW_ARR_F(pbqp_node *, 0);
+ }
+}
+
+static void fill_node_buckets(pbqp *pbqp)
+{
+ unsigned node_index;
+ unsigned node_len;
+
+ assert(pbqp);
+ node_len = pbqp->num_nodes;
+
+ for (node_index = 0; node_index < node_len; ++node_index) {
+ unsigned arity;
+ pbqp_node *node = get_node(pbqp, node_index);
+
+ if (!node) continue;
+
+ arity = ARR_LEN(node->edges);
+
+ /* We have only one bucket for nodes with arity >= 3. */
+ if (arity > 3) {
+ arity = 3;
+ }
+
+ node->bucket_index = ARR_LEN(node_buckets[arity]);
+
+ ARR_APP1(pbqp_node *, node_buckets[arity], node);
+ }
+
+ buckets_filled = 1;
+}
+
+static void normalize_towards_source(pbqp *pbqp, pbqp_edge *edge)
+{
+ pbqp_matrix *mat;
+ pbqp_node *src_node;
+ pbqp_node *tgt_node;
+ vector *src_vec;
+ vector *tgt_vec;
+ int src_len;
+ int tgt_len;
+ int src_index;
+
+ assert(pbqp);
+ assert(edge);
+
+ src_node = edge->src;
+ tgt_node = edge->tgt;
+ assert(src_node);
+ assert(tgt_node);
+
+ src_vec = src_node->costs;
+ tgt_vec = tgt_node->costs;
+ assert(src_vec);
+ assert(tgt_vec);
+
+ src_len = src_vec->len;
+ tgt_len = tgt_vec->len;
+ assert(src_len > 0);
+ assert(tgt_len > 0);
+
+ mat = edge->costs;
+ assert(mat);
+
+ /* Normalize towards source node. */
+ for (src_index = 0; src_index < src_len; ++src_index) {
+ num min = pbqp_matrix_get_row_min(mat, src_index, tgt_vec);
+
+ if (min != 0) {
+ pbqp_matrix_sub_row_value(mat, src_index, tgt_vec, min);
+ src_vec->entries[src_index].data = add(
+ src_vec->entries[src_index].data, min);
+
+ // TODO add to edge_list if inf
+ }
+ }
+}
+
+static void normalize_towards_target(pbqp *pbqp, pbqp_edge *edge)
+{
+ pbqp_matrix *mat;
+ pbqp_node *src_node;
+ pbqp_node *tgt_node;
+ vector *src_vec;
+ vector *tgt_vec;
+ int src_len;
+ int tgt_len;
+ int tgt_index;
+
+ assert(pbqp);
+ assert(edge);
+
+ src_node = edge->src;
+ tgt_node = edge->tgt;
+ assert(src_node);
+ assert(tgt_node);
+
+ src_vec = src_node->costs;
+ tgt_vec = tgt_node->costs;
+ assert(src_vec);
+ assert(tgt_vec);
+
+ src_len = src_vec->len;
+ tgt_len = tgt_vec->len;
+ assert(src_len > 0);
+ assert(tgt_len > 0);
+
+ mat = edge->costs;
+ assert(mat);
+
+ for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
+ num min = pbqp_matrix_get_col_min(mat, tgt_index, src_vec);
+
+ if (min != 0) {
+ pbqp_matrix_sub_col_value(mat, tgt_index, src_vec, min);
+ tgt_vec->entries[tgt_index].data = add(
+ tgt_vec->entries[tgt_index].data, min);
+
+ // TODO add to edge_list if inf
+ }
+ }
+}
+
+static void reorder_node(pbqp_node *node)
+{
+ unsigned arity;
+ unsigned old_arity;
+ unsigned old_bucket_len;
+ unsigned old_bucket_index;
+ pbqp_node **old_bucket;
+ pbqp_node *other;
+
+ if (!buckets_filled) return;
+
+ assert(node);
+
+ arity = ARR_LEN(node->edges);
+
+ /* Same bucket as before */
+ if (arity > 2) return;
+
+ /* Assume node lost one incident edge. */
+ old_arity = arity + 1;
+ old_bucket = node_buckets[old_arity];
+ old_bucket_len = ARR_LEN(old_bucket);
+ old_bucket_index = node->bucket_index;
+
+ if (old_bucket_len <= old_bucket_index ||
+ old_bucket[old_bucket_index] != node) {
+ /* Old arity is new arity, so we have nothing to do. */
+ assert(old_bucket_index < ARR_LEN(node_buckets[arity]) &&
+ node_buckets[arity][old_bucket_index] == node);
+ return;
+ }
+
+ assert(old_bucket[old_bucket_index] == node);
+
+ /* Delete node from old bucket... */
+ other = old_bucket[old_bucket_len - 1];
+ other->bucket_index = old_bucket_index;
+ old_bucket[old_bucket_index] = other;
+ ARR_SHRINKLEN(node_buckets[old_arity], old_bucket_len - 1);
+
+ /* ..and add to new one. */
+ node->bucket_index = ARR_LEN(node_buckets[arity]);
+ ARR_APP1(pbqp_node*, node_buckets[arity], node);
+}
+
+static void simplify_edge(pbqp *pbqp, pbqp_edge *edge)
+{
+ pbqp_matrix *mat;
+ pbqp_node *src_node;
+ pbqp_node *tgt_node;
+ vector *src_vec;
+ vector *tgt_vec;
+ int src_len;
+ int tgt_len;
+
+ assert(pbqp);
+ assert(edge);
+
+ src_node = edge->src;
+ tgt_node = edge->tgt;
+ assert(src_node);
+ assert(tgt_node);
+
+ if (pbqp->dump_file) {
+ char txt[100];
+ sprintf(txt, "Simplification of Edge n%d-n%d", src_node->index, tgt_node->index);
+ dump_section(pbqp->dump_file, 3, txt);
+ }
+
+ src_vec = src_node->costs;
+ tgt_vec = tgt_node->costs;
+ assert(src_vec);
+ assert(tgt_vec);
+
+ src_len = src_vec->len;
+ tgt_len = tgt_vec->len;
+ assert(src_len > 0);
+ assert(tgt_len > 0);
+
+ mat = edge->costs;
+ assert(mat);
+
+ if (pbqp->dump_file) {
+ fputs("Input:<br>\n", pbqp->dump_file);
+ dump_simplifyedge(pbqp, edge);
+ }
+
+ normalize_towards_source(pbqp, edge);
+ normalize_towards_target(pbqp, edge);
+
+ if (pbqp->dump_file) {
+ fputs("<br>\nOutput:<br>\n", pbqp->dump_file);
+ dump_simplifyedge(pbqp, edge);
+ }
+
+ if (pbqp_matrix_is_zero(mat, src_vec, tgt_vec)) {
+ if (pbqp->dump_file) {
+ fputs("edge has been eliminated", pbqp->dump_file);
+
+ delete_edge(edge);
+ reorder_node(src_node);
+ reorder_node(tgt_node);
+ }
+ }
+}