Add functions to find minimum cut from a maximum flow

src/flow.c

@@ -1,7 +1,6 @@
 #include <malloc.h>
 #include <float.h>
 #include "gtsp.h"
-#include "flow.h"
 #include "util.h"
 
 int flow_find_max_flow(
@@ -58,6 +57,9 @@ int flow_find_max_flow(
         }
     }
 
+    rval = flow_mark_reachable_nodes(digraph, residual_caps, from);
+    abort_if(rval, "flow_mark_reachable_nodes failed");
+
     CLEANUP:
     if (path_edges) free(path_edges);
     if (residual_caps) free(residual_caps);
@@ -151,4 +153,44 @@ int flow_find_augmenting_path(
     if (parent_edges) free(parent_edges);
     if (queue) free(queue);
     return rval;
+}
+
+int flow_mark_reachable_nodes(
+        struct Graph *graph, double *residual_caps, struct Node *from)
+{
+    int rval = 0;
+
+    struct Node **stack;
+    int stack_top = 0;
+
+    stack = (struct Node**) malloc(graph->node_count * sizeof(struct Node*));
+    abort_if(!stack, "could not allocate stack");
+
+    for (int i = 0; i < graph->node_count; i++)
+        graph->nodes[i].mark = 0;
+
+    from->mark = 1;
+    stack[stack_top++] = from;
+
+    while(stack_top > 0)
+    {
+        struct Node *n = stack[--stack_top];
+
+        for (int j = 0; j < n->degree; j++)
+        {
+            struct Edge *e = n->adj[j].edge;
+            struct Node *neighbor = n->adj[j].neighbor;
+
+            if(neighbor->mark) continue;
+            if(residual_caps[e->index] <= 0) continue;
+
+            stack[stack_top++] = neighbor;
+            neighbor->mark = 1;
+        }
+
+    }
+
+    CLEANUP:
+    if(stack) free(stack);
+    return rval;
 }

src/graph.c

@@ -52,7 +52,7 @@ int graph_build(
         int a = edges[2 * i];
         int b = edges[2 * i + 1];
         graph->nodes[a].degree++;
-        if(!is_directed) graph->nodes[b].degree++;
+        if (!is_directed) graph->nodes[b].degree++;
 
         graph->edges[i].reverse = 0;
         graph->edges[i].index = i;
@@ -80,7 +80,7 @@ int graph_build(
         n->adj[n->degree].edge = &graph->edges[i];
         n->degree++;
 
-        if(!is_directed)
+        if (!is_directed)
         {
             n = &graph->nodes[b];
             n->adj[n->degree].neighbor_index = a;
@@ -179,3 +179,20 @@ void get_delta(
     for (int i = 0; i < island_node_count; i++)
         marks[island_nodes[i]] = 0;
 }
+
+int get_cut_edges_from_marks(
+        struct Graph *graph, int *cut_edges_count, struct Edge **cut_edges)
+{
+    *cut_edges_count = 0;
+
+    for (int i = 0; i < graph->edge_count; ++i)
+    {
+        struct Edge *e = &graph->edges[i];
+        struct Node *from = e->from;
+        struct Node *to = e->to;
+        if (from->mark && !to->mark)
+            cut_edges[(*cut_edges_count)++] = e;
+    }
+
+    return 0;
+}

src/graph.h

@@ -70,4 +70,7 @@ void get_delta(
 int graph_build_directed_from_undirected
         (const struct Graph *graph, struct Graph *digraph);
 
+int get_cut_edges_from_marks(
+        struct Graph *graph, int *cut_edges_count, struct Edge **cut_edges);
+
 #endif

src/main.c

@@ -29,6 +29,8 @@ int test_max_flow()
     double *flow = 0;
     double flow_value;
 
+    struct Edge **cut_edges = 0;
+
     FILE *f = fopen("tmp/flow.in", "r");
     abort_if(!f, "could not open input file");
 
@@ -53,8 +55,8 @@ int test_max_flow()
         rval = fscanf(f, "%d %d %d ", &from, &to, &cap);
         abort_if(rval != 3, "invalid input format");
 
-        edges[i*4] = edges[i*4+3] = from;
-        edges[i*4+1] = edges[i*4 + 2] = to;
+        edges[i * 4] = edges[i * 4 + 3] = from;
+        edges[i * 4 + 1] = edges[i * 4 + 2] = to;
         capacities[2 * i] = cap;
         capacities[2 * i + 1] = 0;
     }
@@ -64,8 +66,8 @@ int test_max_flow()
 
     for (int i = 0; i < edge_count; i++)
     {
-        graph.edges[2*i].reverse = &graph.edges[2*i+1];
-        graph.edges[2*i+1].reverse = &graph.edges[2*i];
+        graph.edges[2 * i].reverse = &graph.edges[2 * i + 1];
+        graph.edges[2 * i + 1].reverse = &graph.edges[2 * i];
     }
 
     flow = (double *) malloc(graph.edge_count * sizeof(double));
@@ -81,12 +83,38 @@ int test_max_flow()
     for (int i = 0; i < graph.edge_count; i++)
     {
         struct Edge *e = &graph.edges[i];
-        if(flow[e->index] <= 0) continue;
+        if (flow[e->index] <= 0) continue;
+
+        log_info("  %d %d %6.2f / %6.2f\n", e->from->index, e->to->index,
+                flow[e->index], capacities[e->index]);
+    }
+
+    log_info("Nodes reachable from origin on residual graph:\n");
+    for (int i = 0; i < graph.node_count; i++)
+    {
+        struct Node *n = &graph.nodes[i];
+        if(n->mark)
+            log_info("  %d\n", n->index);
+    }
+
+    int cut_edges_count = 0;
+    cut_edges = (struct Edge**) malloc(graph.edge_count * sizeof(struct Edge*));
+    abort_if(!cut_edges, "could not allocate cut_edges");
+
+    rval = get_cut_edges_from_marks(&graph, &cut_edges_count, cut_edges);
+    abort_if(rval, "get_cut_edges_from_marks failed");
 
-        log_info("  %d %d %6.2f / %6.2f\n", e->from->index, e->to->index, flow[e->index], capacities[e->index]);
+    log_info("Min cut edges:\n");
+    for (int i = 0; i < cut_edges_count; i++)
+    {
+        struct Edge *e = cut_edges[i];
+        if(capacities[e->index] <= 0) continue;
+        log_info("  %d %d\n", e->from->index, e->to->index);
     }
 
+
     CLEANUP:
+    if(cut_edges) free(cut_edges);
     if (capacities) free(capacities);
     if (edges) free(edges);
     if (flow) free(flow);