askrene: add a simple MCF solver

Changelog-EXPERIMENTAL: askrene: add a simple MCF solver

Signed-off-by: Lagrang3 <lagrang3@protonmail.com>

plugins/askrene/algorithm.c

@@ -322,3 +322,87 @@ s64 node_balance(const struct graph *graph,
 	}
 	return balance;
 }
+
+
+bool simple_mcf(const tal_t *ctx, const struct graph *graph,
+		const struct node source, const struct node destination,
+		s64 *capacity, s64 amount, const s64 *cost)
+{
+	tal_t *this_ctx = tal(ctx, tal_t);
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+	const size_t max_num_nodes = graph_max_num_nodes(graph);
+	s64 remaining_amount = amount;
+
+	if (amount < 0)
+		goto finish;
+
+	if (!graph || source.idx >= max_num_nodes ||
+	    destination.idx >= max_num_nodes || !capacity || !cost)
+		goto finish;
+
+	if (tal_count(capacity) != max_num_arcs ||
+	    tal_count(cost) != max_num_arcs)
+		goto finish;
+
+	struct arc *prev = tal_arr(this_ctx, struct arc, max_num_nodes);
+	s64 *distance = tal_arrz(this_ctx, s64, max_num_nodes);
+	s64 *potential = tal_arrz(this_ctx, s64, max_num_nodes);
+
+	if (!prev || !distance || !potential)
+		goto finish;
+
+	/* FIXME: implement this algorithm as a search for matching negative and
+	 * positive balance nodes, so that we can use it to adapt a flow
+	 * structure for changes in the cost function. */
+	while (remaining_amount > 0) {
+		if (!dijkstra_path(this_ctx, graph, source, destination,
+				   /* prune = */ true, capacity, 1, cost,
+				   potential, prev, distance))
+			goto finish;
+
+		/* traverse the path and see how much flow we can send */
+		s64 delta = get_augmenting_flow(graph, source, destination,
+						capacity, prev);
+
+		/* commit that flow to the path */
+		delta = MIN(remaining_amount, delta);
+		assert(delta > 0 && delta <= remaining_amount);
+
+		augment_flow(graph, source, destination, prev, capacity, delta);
+		remaining_amount -= delta;
+
+		/* update potentials */
+		for (u32 n = 0; n < max_num_nodes; n++) {
+			/* see page 323 of Ahuja-Magnanti-Orlin.
+			 * Whether we prune or not the Dijkstra search, the
+			 * following potentials will keep reduced costs
+			 * non-negative. */
+			potential[n] -=
+			    MIN(distance[destination.idx], distance[n]);
+		}
+	}
+finish:
+	tal_free(this_ctx);
+	return remaining_amount == 0;
+}
+
+s64 flow_cost(const struct graph *graph, const s64 *capacity, const s64 *cost)
+{
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+	s64 total_cost = 0;
+
+	assert(graph && capacity && cost);
+	assert(tal_count(capacity) == max_num_arcs &&
+	       tal_count(cost) == max_num_arcs);
+
+	for (u32 i = 0; i < max_num_arcs; i++) {
+		struct arc arc = {.idx = i};
+		struct arc dual = arc_dual(graph, arc);
+
+		if (arc_is_dual(graph, arc))
+			continue;
+
+		total_cost += capacity[dual.idx] * cost[arc.idx];
+	}
+	return total_cost;
+}

plugins/askrene/algorithm.h

@@ -115,4 +115,42 @@ bool simple_feasibleflow(const tal_t *ctx, const struct graph *graph,
 s64 node_balance(const struct graph *graph, const struct node node,
 		 const s64 *capacity);
 
+
+/* Finds the minimum cost flow that satisfy the capacity constraints:
+ * 	flow[i] <= capacity[i]
+ * and supply/demand constraints:
+ * 	supply[source] = demand[destination] = amount
+ * 	supply/demand[node] = 0 for every other node
+ *
+ * It uses successive shortest path algorithm.
+ *
+ * input:
+ * @ctx: tal context for internal allocation
+ * @graph: topological information of the graph
+ * @source: source node
+ * @destination: destination node
+ * @capacity: arcs capacity
+ * @amount: desired balance at the destination
+ * @cost: cost per unit of flow
+ *
+ * output:
+ * @capacity: residual capacity
+ * returns true if the balance constraint can be satisfied
+ *
+ * precondition:
+ * |capacity|=graph_max_num_arcs
+ * |cost|=graph_max_num_arcs
+ * amount>=0
+ * */
+bool simple_mcf(const tal_t *ctx, const struct graph *graph,
+		const struct node source, const struct node destination,
+		s64 *capacity, s64 amount, const s64 *cost);
+
+/* Compute the cost of a flow in the network.
+ *
+ * @graph: network topology
+ * @capacity: residual capacity (encodes the flow)
+ * @cost: cost per unit of flow */
+s64 flow_cost(const struct graph *graph, const s64 *capacity, const s64 *cost);
+
 #endif /* LIGHTNING_PLUGINS_ASKRENE_ALGORITHM_H */

plugins/askrene/test/Makefile

@@ -10,7 +10,7 @@ $(PLUGIN_RENEPAY_TEST_OBJS): $(PLUGIN_ASKRENE_SRC)
 
 PLUGIN_ASKRENE_TEST_COMMON_OBJS :=
 
-plugins/askrene/test/run-bfs plugins/askrene/test/run-dijkstra plugins/askrene/test/run-flow: \
+plugins/askrene/test/run-bfs plugins/askrene/test/run-dijkstra plugins/askrene/test/run-flow plugins/askrene/test/run-mcf: \
 	plugins/askrene/priorityqueue.o \
 	plugins/askrene/graph.o
 

plugins/askrene/test/run-mcf.c

@@ -0,0 +1,68 @@
+#include "config.h"
+#include <assert.h>
+#include <ccan/tal/tal.h>
+#include <common/setup.h>
+#include <inttypes.h>
+#include <plugins/askrene/graph.h>
+#include <stdio.h>
+
+#include "../algorithm.c"
+
+#define CHECK(arg) if(!(arg)){fprintf(stderr, "failed CHECK at line %d: %s\n", __LINE__, #arg); abort();}
+
+#define MAX_NODES 256
+#define MAX_ARCS 256
+#define DUAL_BIT 7
+
+int main(int argc, char *argv[])
+{
+	common_setup(argv[0]);
+	printf("Allocating a memory context\n");
+	tal_t *ctx = tal(NULL, tal_t);
+	assert(ctx);
+
+	printf("Allocating a graph\n");
+	struct graph *graph = graph_new(ctx, MAX_NODES, MAX_ARCS, DUAL_BIT);
+	assert(graph);
+
+	s64 *capacity = tal_arrz(ctx, s64, MAX_ARCS);
+	s64 *cost = tal_arrz(ctx, s64, MAX_ARCS);
+
+	graph_add_arc(graph, arc_obj(0), node_obj(0), node_obj(1));
+	capacity[0] = 2, cost[0] = 0;
+	graph_add_arc(graph, arc_obj(1), node_obj(0), node_obj(2));
+	capacity[1] = 2, cost[1] = 0;
+	graph_add_arc(graph, arc_obj(2), node_obj(1), node_obj(3));
+	capacity[2] = 1, cost[2] = 1;
+	graph_add_arc(graph, arc_obj(3), node_obj(1), node_obj(4));
+	capacity[3] = 1, cost[3] = 2;
+	graph_add_arc(graph, arc_obj(4), node_obj(2), node_obj(3));
+	capacity[4] = 2, cost[4] = 1;
+	graph_add_arc(graph, arc_obj(5), node_obj(2), node_obj(4));
+	capacity[5] = 1, cost[5] = 2;
+	graph_add_arc(graph, arc_obj(6), node_obj(3), node_obj(5));
+	capacity[6] = 3, cost[6] = 0;
+	graph_add_arc(graph, arc_obj(7), node_obj(4), node_obj(5));
+	capacity[7] = 3, cost[7] = 0;
+
+	struct node src = {.idx = 0};
+	struct node dst = {.idx = 5};
+
+	bool result = simple_mcf(ctx, graph, src, dst, capacity, 4, cost);
+	CHECK(result);
+
+	CHECK(node_balance(graph, src, capacity) == -4);
+	CHECK(node_balance(graph, dst, capacity) == 4);
+
+	for (u32 i = 1; i < 4; i++)
+		CHECK(node_balance(graph, node_obj(i), capacity) == 0);
+
+	const s64 total_cost = flow_cost(graph, capacity, cost);
+	printf("best flow cost: %" PRIi64 "\n", total_cost);
+	CHECK(total_cost == 5);
+
+	printf("Freeing memory\n");
+	ctx = tal_free(ctx);
+	common_shutdown();
+	return 0;
+}