askrene: move routines only accessed by the child process into child/.

diff --git a/plugins/askrene/Makefile b/plugins/askrene/Makefile
index 9e210f08..3c3f5c4b 100644
--- a/plugins/askrene/Makefile
+++ b/plugins/askrene/Makefile
@@ -1,32 +1,22 @@
-PLUGIN_ASKRENE_SRC := \
+PLUGIN_ASKRENE_PARENT_SRC := \
 	plugins/askrene/askrene.c \
 	plugins/askrene/datastore_wire.c \
 	plugins/askrene/layer.c \
 	plugins/askrene/reserve.c \
-	plugins/askrene/mcf.c \
-	plugins/askrene/dijkstra.c \
-	plugins/askrene/flow.c \
-	plugins/askrene/refine.c \
-	plugins/askrene/explain_failure.c \
-	plugins/askrene/graph.c \
-	plugins/askrene/priorityqueue.c \
-	plugins/askrene/algorithm.c \
+
+PLUGIN_ASKRENE_CHILD_SRC := \
+	plugins/askrene/child/mcf.c \
+	plugins/askrene/child/dijkstra.c \
+	plugins/askrene/child/flow.c \
+	plugins/askrene/child/refine.c \
+	plugins/askrene/child/explain_failure.c \
+	plugins/askrene/child/graph.c \
+	plugins/askrene/child/priorityqueue.c \
+	plugins/askrene/child/algorithm.c \
 	plugins/askrene/child/child_log.c \
 
-PLUGIN_ASKRENE_HEADER := \
-	plugins/askrene/askrene.h \
-	plugins/askrene/datastore_wire.h \
-	plugins/askrene/layer.h \
-	plugins/askrene/reserve.h \
-	plugins/askrene/mcf.h \
-	plugins/askrene/dijkstra.h \
-	plugins/askrene/flow.h \
-	plugins/askrene/refine.h \
-	plugins/askrene/explain_failure.h \
-	plugins/askrene/graph.h \
-	plugins/askrene/priorityqueue.h \
-	plugins/askrene/algorithm.h \
-	plugins/askrene/child/child_log.h \
+PLUGIN_ASKRENE_SRC := $(PLUGIN_ASKRENE_PARENT_SRC) $(PLUGIN_ASKRENE_CHILD_SRC)
+PLUGIN_ASKRENE_HEADER := $(PLUGIN_ASKRENE_SRC:.c=.h)
 
 PLUGIN_ASKRENE_OBJS := $(PLUGIN_ASKRENE_SRC:.c=.o)
 
diff --git a/plugins/askrene/algorithm.c b/plugins/askrene/algorithm.c
deleted file mode 100644
index d253325a..00000000
--- a/plugins/askrene/algorithm.c
+++ /dev/null
@@ -1,670 +0,0 @@
-#include "config.h"
-#include <ccan/bitmap/bitmap.h>
-#include <ccan/tal/tal.h>
-#include <plugins/askrene/algorithm.h>
-#include <plugins/askrene/priorityqueue.h>
-
-static const s64 INFINITE = INT64_MAX;
-
-#define MAX(x, y) (((x) > (y)) ? (x) : (y))
-#define MIN(x, y) (((x) < (y)) ? (x) : (y))
-
-bool BFS_path(const tal_t *ctx, const struct graph *graph,
-	      const struct node source, const struct node destination,
-	      const s64 *capacity, const s64 cap_threshold, struct arc *prev)
-{
-	const tal_t *this_ctx = tal(ctx, tal_t);
-	bool target_found = false;
-	assert(graph);
-	const size_t max_num_arcs = graph_max_num_arcs(graph);
-	const size_t max_num_nodes = graph_max_num_nodes(graph);
-
-	/* check preconditions */
-	assert(source.idx < max_num_nodes);
-	assert(capacity);
-	assert(prev);
-	assert(tal_count(capacity) == max_num_arcs);
-	assert(tal_count(prev) == max_num_nodes);
-
-	for (size_t i = 0; i < max_num_nodes; i++)
-		prev[i].idx = INVALID_INDEX;
-
-	/* A minimalistic queue is implemented here. Nodes are not visited more
-	 * than once, therefore a maximum size of max_num_nodes is sufficient.
-	 * max_num_arcs would work as well but we expect max_num_arcs to be a
-	 * factor >10 greater than max_num_nodes. */
-	u32 *queue = tal_arr(this_ctx, u32, max_num_nodes);
-	size_t queue_start = 0, queue_end = 0;
-
-	queue[queue_end++] = source.idx;
-
-	while (queue_start < queue_end) {
-		struct node cur = {.idx = queue[queue_start++]};
-
-		if (cur.idx == destination.idx) {
-			target_found = true;
-			break;
-		}
-
-		for (struct arc arc = node_adjacency_begin(graph, cur);
-		     !node_adjacency_end(arc);
-		     arc = node_adjacency_next(graph, arc)) {
-			/* check if this arc is traversable */
-			if (capacity[arc.idx] < cap_threshold)
-				continue;
-
-			const struct node next = arc_head(graph, arc);
-
-			/* if that node has been seen previously */
-			if (prev[next.idx].idx != INVALID_INDEX ||
-			    next.idx == source.idx)
-				continue;
-
-			prev[next.idx] = arc;
-
-			assert(queue_end < max_num_nodes);
-			queue[queue_end++] = next.idx;
-		}
-	}
-
-	tal_free(this_ctx);
-	return target_found;
-}
-
-bool dijkstra_path(const tal_t *ctx, const struct graph *graph,
-		   const struct node source, const struct node destination,
-		   bool prune, const s64 *capacity, const s64 cap_threshold,
-		   const s64 *cost, const s64 *potential, struct arc *prev,
-		   s64 *distance)
-{
-	bool target_found = false;
-	assert(graph);
-	const size_t max_num_arcs = graph_max_num_arcs(graph);
-	const size_t max_num_nodes = graph_max_num_nodes(graph);
-	const tal_t *this_ctx = tal(ctx, tal_t);
-
-	/* check preconditions */
-	assert(source.idx<max_num_nodes);
-	assert(cost);
-	assert(capacity);
-	assert(prev);
-	assert(distance);
-
-	/* if prune is true then the destination cannot be invalid */
-	assert(destination.idx < max_num_nodes || !prune);
-
-	assert(tal_count(cost) == max_num_arcs);
-	assert(tal_count(capacity) == max_num_arcs);
-	assert(tal_count(prev) == max_num_nodes);
-	assert(tal_count(distance) == max_num_nodes);
-
-	/* FIXME: maybe this is unnecessary */
-	bitmap *visited = tal_arrz(this_ctx, bitmap,
-				   BITMAP_NWORDS(max_num_nodes));
-
-	for (size_t i = 0; i < max_num_nodes; ++i)
-		prev[i].idx = INVALID_INDEX;
-
-	struct priorityqueue *q;
-	q = priorityqueue_new(this_ctx, max_num_nodes);
-	const s64 *const dijkstra_distance = priorityqueue_value(q);
-
-	priorityqueue_init(q);
-	priorityqueue_update(q, source.idx, 0);
-
-	while (!priorityqueue_empty(q)) {
-		const u32 cur = priorityqueue_top(q);
-		priorityqueue_pop(q);
-
-		/* FIXME: maybe this is unnecessary */
-		if (bitmap_test_bit(visited, cur))
-			continue;
-		bitmap_set_bit(visited, cur);
-
-		if (cur == destination.idx) {
-			target_found = true;
-			if (prune)
-				break;
-		}
-
-		for (struct arc arc =
-			 node_adjacency_begin(graph, node_obj(cur));
-		     !node_adjacency_end(arc);
-		     arc = node_adjacency_next(graph, arc)) {
-			/* check if this arc is traversable */
-			if (capacity[arc.idx] < cap_threshold)
-				continue;
-
-			const struct node next = arc_head(graph, arc);
-
-			const s64 cij = cost[arc.idx] - potential[cur] +
-					potential[next.idx];
-
-			/* Dijkstra only works with non-negative weights */
-			assert(cij >= 0);
-
-			if (dijkstra_distance[next.idx] <=
-			    dijkstra_distance[cur] + cij)
-				continue;
-
-			priorityqueue_update(q, next.idx,
-					     dijkstra_distance[cur] + cij);
-			prev[next.idx] = arc;
-		}
-	}
-	for (size_t i = 0; i < max_num_nodes; i++)
-		distance[i] = dijkstra_distance[i];
-
-	tal_free(this_ctx);
-	return target_found;
-}
-
-/* Get the max amount of flow one can send from source to target along the path
- * encoded in `prev`. */
-static s64 get_augmenting_flow(const struct graph *graph,
-			       const struct node source,
-			       const struct node target, const s64 *capacity,
-			       const struct arc *prev)
-{
-	const size_t max_num_nodes = graph_max_num_nodes(graph);
-	const size_t max_num_arcs = graph_max_num_arcs(graph);
-	assert(max_num_nodes == tal_count(prev));
-	assert(max_num_arcs == tal_count(capacity));
-
-	/* count the number of arcs in the path */
-	int path_length = 0;
-	s64 flow = INFINITE;
-
-	struct node cur = target;
-	while (cur.idx != source.idx) {
-		assert(cur.idx < max_num_nodes);
-		const struct arc arc = prev[cur.idx];
-		assert(arc.idx < max_num_arcs);
-		flow = MIN(flow, capacity[arc.idx]);
-
-		/* we are traversing in the opposite direction to the flow,
-		 * hence the next node is at the tail of the arc. */
-		cur = arc_tail(graph, arc);
-
-		/* We may never have a path exceeds the number of nodes, it this
-		 * happens it means we have an infinite loop. */
-		path_length++;
-		if(path_length >= max_num_nodes){
-			flow = -1;
-			break;
-		}
-	}
-
-	assert(flow < INFINITE && flow > 0);
-	return flow;
-}
-
-
-/* Helper.
- * Sends an amount of flow through an arc, changing the flow balance of the
- * nodes connected by the arc and the [residual] capacity of the arc and its
- * dual. */
-static void sendflow(const struct graph *graph, const struct arc arc,
-		     const s64 flow, s64 *arc_capacity, s64 *node_balance)
-{
-	const struct arc dual = arc_dual(graph, arc);
-
-	arc_capacity[arc.idx] -= flow;
-	arc_capacity[dual.idx] += flow;
-
-	if (node_balance) {
-		const struct node src = arc_tail(graph, arc),
-				  dst = arc_tail(graph, dual);
-
-		node_balance[src.idx] -= flow;
-		node_balance[dst.idx] += flow;
-	}
-}
-
-/* Augment a `flow` amount along the path defined by `prev`.*/
-static void augment_flow(const struct graph *graph,
-			 const struct node source,
-			 const struct node target,
-			 const struct arc *prev,
-			 s64 *excess,
-			 s64 *capacity,
-			 s64 flow)
-{
-	const size_t max_num_nodes = graph_max_num_nodes(graph);
-	const size_t max_num_arcs = graph_max_num_arcs(graph);
-	assert(max_num_nodes == tal_count(prev));
-	assert(max_num_arcs == tal_count(capacity));
-
-	struct node cur = target;
-	/* count the number of arcs in the path */
-	int path_length = 0;
-
-	while (cur.idx != source.idx) {
-		assert(cur.idx < max_num_nodes);
-		const struct arc arc = prev[cur.idx];
-
-		sendflow(graph, arc, flow, capacity, excess);
-
-		/* we are traversing in the opposite direction to the flow,
-		 * hence the next node is at the tail of the arc. */
-		cur = arc_tail(graph, arc);
-
-		/* We may never have a path exceeds the number of nodes, it this
-		 * happens it means we have an infinite loop. */
-		path_length++;
-		if (path_length >= max_num_nodes)
-			break;
-	}
-	assert(path_length < max_num_nodes);
-}
-
-bool simple_feasibleflow(const tal_t *ctx,
-			 const struct graph *graph,
-			 const struct node source,
-			 const struct node destination,
-			 s64 *capacity,
-			 s64 amount)
-{
-	const tal_t *this_ctx = tal(ctx, tal_t);
-	assert(graph);
-	const size_t max_num_arcs = graph_max_num_arcs(graph);
-	const size_t max_num_nodes = graph_max_num_nodes(graph);
-
-	/* check preconditions */
-	assert(amount > 0);
-	assert(source.idx < max_num_nodes);
-	assert(destination.idx < max_num_nodes);
-	assert(capacity);
-	assert(tal_count(capacity) == max_num_arcs);
-
-	/* path information
-	 * prev: is the id of the arc that lead to the node. */
-	struct arc *prev = tal_arr(this_ctx, struct arc, max_num_nodes);
-	if (!prev)
-		goto finish;
-
-	while (amount > 0) {
-		/* find a path from source to target */
-		if (!BFS_path(this_ctx, graph, source, destination, capacity, 1,
-			      prev))
-			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(amount, delta);
-		assert(delta > 0 && delta <= amount);
-
-		augment_flow(graph, source, destination, prev, NULL, capacity,
-			     delta);
-		amount -= delta;
-	}
-finish:
-	tal_free(this_ctx);
-	return amount == 0;
-}
-
-s64 node_balance(const struct graph *graph,
-		 const struct node node,
-		 const s64 *capacity)
-{
-	s64 balance = 0;
-
-	for (struct arc arc = node_adjacency_begin(graph, node);
-	     !node_adjacency_end(arc); arc = node_adjacency_next(graph, arc)) {
-		struct arc dual = arc_dual(graph, arc);
-
-		if (arc_is_dual(graph, arc))
-			balance += capacity[arc.idx];
-		else
-			balance -= capacity[dual.idx];
-	}
-	return balance;
-}
-
-/* Helper.
- * Compute the reduced cost of an arc. */
-static s64 reduced_cost(const struct graph *graph, const struct arc arc,
-			const s64 *cost, const s64 *potential)
-{
-	struct node src = arc_tail(graph, arc);
-	struct node dst = arc_head(graph, arc);
-	return cost[arc.idx] - potential[src.idx] + potential[dst.idx];
-}
-
-/* Finds an optimal path from the source to the nearest sink node, by definition
- * a node i is a sink if node_balance[i]<0. It uses a reduced cost:
- *	reduced_cost[i,j] = cost[i,j] - potential[i] + potential[j]
- *
- * */
-static struct node dijkstra_nearest_sink(const tal_t *ctx,
-					 const struct graph *graph,
-					 const struct node source,
-					 const s64 *node_balance,
-					 const s64 *capacity,
-					 const s64 cap_threshold,
-					 const s64 *cost,
-					 const s64 *potential,
-					 struct arc *prev,
-					 s64 *distance)
-{
-	struct node target = {.idx = INVALID_INDEX};
-	const tal_t *this_ctx = tal(ctx, tal_t);
-
-	/* check preconditions */
-	assert(graph);
-	assert(node_balance);
-	assert(capacity);
-	assert(cost);
-	assert(potential);
-	assert(prev);
-	assert(distance);
-
-	const size_t max_num_arcs = graph_max_num_arcs(graph);
-	const size_t max_num_nodes = graph_max_num_nodes(graph);
-
-	assert(source.idx < max_num_nodes);
-	assert(tal_count(node_balance) == max_num_nodes);
-	assert(tal_count(capacity) == max_num_arcs);
-	assert(tal_count(cost) == max_num_arcs);
-	assert(tal_count(potential) == max_num_nodes);
-	assert(tal_count(prev) == max_num_nodes);
-	assert(tal_count(distance) == max_num_nodes);
-
-	for (size_t i = 0; i < max_num_arcs; i++) {
-		/* is this arc saturated? */
-		if (capacity[i] < cap_threshold)
-			continue;
-
-		struct arc arc = {.idx = i};
-		struct node tail = arc_tail(graph, arc);
-		struct node head = arc_head(graph, arc);
-		s64 red_cost =
-		    cost[i] - potential[tail.idx] + potential[head.idx];
-
-		/* reducted cost cannot be negative for non saturated arcs,
-		 * otherwise Dijkstra does not work. */
-		if (red_cost < 0)
-			goto finish;
-	}
-
-	for (size_t i = 0; i < max_num_nodes; ++i)
-		prev[i].idx = INVALID_INDEX;
-
-/* Only in debug mode we keep track of visited nodes. */
-#ifdef ASKRENE_UNITTEST
-	bitmap *visited =
-	    tal_arrz(this_ctx, bitmap, BITMAP_NWORDS(max_num_nodes));
-#endif
-
-	struct priorityqueue *q;
-	q = priorityqueue_new(this_ctx, max_num_nodes);
-	const s64 *const dijkstra_distance = priorityqueue_value(q);
-
-	priorityqueue_init(q);
-	priorityqueue_update(q, source.idx, 0);
-
-	while (!priorityqueue_empty(q)) {
-		const u32 idx = priorityqueue_top(q);
-		const struct node cur = {.idx = idx};
-		priorityqueue_pop(q);
-
-/* Only in debug mode we keep track of visited nodes. */
-#ifdef ASKRENE_UNITTEST
-		assert(!bitmap_test_bit(visited, cur.idx));
-		bitmap_set_bit(visited, cur.idx);
-#endif
-
-		if (node_balance[cur.idx] < 0) {
-			target = cur;
-			break;
-		}
-
-		for (struct arc arc = node_adjacency_begin(graph, cur);
-		     !node_adjacency_end(arc);
-		     arc = node_adjacency_next(graph, arc)) {
-			/* check if this arc is traversable */
-			if (capacity[arc.idx] < cap_threshold)
-				continue;
-
-			const struct node next = arc_head(graph, arc);
-
-			const s64 cij = cost[arc.idx] - potential[cur.idx] +
-					potential[next.idx];
-
-			/* Dijkstra only works with non-negative weights */
-			assert(cij >= 0);
-
-			if (dijkstra_distance[next.idx] <=
-			    dijkstra_distance[cur.idx] + cij)
-				continue;
-
-			priorityqueue_update(q, next.idx,
-					     dijkstra_distance[cur.idx] + cij);
-			prev[next.idx] = arc;
-		}
-	}
-	for (size_t i = 0; i < max_num_nodes; i++)
-		distance[i] = dijkstra_distance[i];
-
-finish:
-	tal_free(this_ctx);
-	return target;
-}
-
-/* Problem: find a potential and capacity redistribution such that:
- *	excess[all nodes] = 0
- *	capacity[all arcs] >= 0
- *	cost/potential [i,j] < 0 implies capacity[i,j] = 0
- *
- *	Q. Is this a feasible solution?
- *
- *	A. If we use flow conserving function sendflow, then
- *	if for all nodes excess[i] = 0 and capacity[i,j] >= 0 for all arcs
- *	then we have reached a feasible flow.
- *
- *	Q. Is this flow optimal?
- *
- *	A. According to Theorem 9.4 (Ahuja page 309) we have reached an optimal
- *	solution if we are able to find a potential and flow that satisfy the
- *	slackness optimality conditions:
- *
- *		if cost_reduced[i,j] > 0 then x[i,j] = 0
- *		if 0 < x[i,j] < u[i,j] then cost_reduced[i,j] = 0
- *		if cost_reduced[i,j] < 0 then x[i,j] = u[i,j]
- *
- *	In our representation the slackness optimality conditions are equivalent
- *	to the following condition in the residual network:
- *
- *		cost_reduced[i,j] < 0 then capacity[i,j] = 0
- *
- *	Therefore yes, the solution is optimal.
- *
- *	Q. Why is this useful?
- *
- *	A. It can be used to compute a MCF from scratch or build an optimal
- *	solution starting from a non-optimal one, eg. if we first test the
- *	solution feasibility we already have a solution canditate, we use that
- *	flow as input to this function, in another example we might have an
- *	algorithm that changes the cost function at every iteration and we need
- *	to find the MCF every time.
- * */
-bool mcf_refinement(const tal_t *ctx,
-		    const struct graph *graph,
-		    s64 *excess,
-		    s64 *capacity,
-		    const s64 *cost,
-		    s64 *potential)
-{
-	bool solved = false;
-	const tal_t *this_ctx = tal(ctx, tal_t);
-
-	assert(graph);
-	assert(excess);
-	assert(capacity);
-	assert(cost);
-	assert(potential);
-
-	const size_t max_num_arcs = graph_max_num_arcs(graph);
-	const size_t max_num_nodes = graph_max_num_nodes(graph);
-
-	assert(tal_count(excess) == max_num_nodes);
-	assert(tal_count(capacity) == max_num_arcs);
-	assert(tal_count(cost) == max_num_arcs);
-	assert(tal_count(potential) == max_num_nodes);
-
-	s64 total_excess = 0;
-	for (u32 i = 0; i < max_num_nodes; i++)
-		total_excess += excess[i];
-
-	if (total_excess)
-		/* there is no way to satisfy the constraints if supply does not
-		 * match demand */
-		goto finish;
-
-	/* Enforce the complementary slackness condition, rolls back
-	 * constraints.  */
-	for (u32 arc_id = 0; arc_id < max_num_arcs; arc_id++) {
-		struct arc arc = {.idx = arc_id};
-		if(!arc_enabled(graph, arc))
-			continue;
-		const s64 r = capacity[arc.idx];
-		if (reduced_cost(graph, arc, cost, potential) < 0 && r > 0) {
-			/* This arc's reduced cost is negative and non
-			 * saturated. */
-			sendflow(graph, arc, r, capacity, excess);
-		}
-	}
-
-	struct arc *prev = tal_arr(this_ctx, struct arc, max_num_nodes);
-	s64 *distance = tal_arrz(this_ctx, s64, max_num_nodes);
-	if (!prev || !distance)
-		goto finish;
-
-	/* Now build back constraints again keeping the complementary slackness
-	 * condition. */
-	for (u32 node_id = 0; node_id < max_num_nodes; node_id++) {
-		struct node src = {.idx = node_id};
-
-		/* is this node a source */
-		while (excess[src.idx] > 0) {
-
-			/* where is the nearest sink */
-			struct node dst = dijkstra_nearest_sink(
-			    this_ctx, graph, src, excess, capacity, 1, cost,
-			    potential, prev, distance);
-
-			if (dst.idx >= max_num_nodes)
-				/* we failed to find a reacheable sink */
-				goto finish;
-
-			/* traverse the path and see how much flow we can send
-			 */
-			s64 delta = get_augmenting_flow(graph, src, dst,
-							capacity, prev);
-
-			delta = MIN(excess[src.idx], delta);
-			delta = MIN(-excess[dst.idx], delta);
-			assert(delta > 0);
-
-			/* commit that flow to the path */
-			augment_flow(graph, src, dst, prev, excess, capacity,
-				     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[dst.idx], distance[n]);
-			}
-		}
-	}
-
-#ifdef ASKRENE_UNITTEST
-	/* verify that we have satisfied all constraints */
-	for (u32 i = 0; i < max_num_nodes; i++) {
-		assert(excess[i] == 0);
-	}
-	for (u32 i = 0; i < max_num_arcs; i++) {
-		struct arc arc = {.idx = i};
-		if(!arc_enabled(graph, arc))
-			continue;
-		const s64 cap = capacity[arc.idx];
-		const s64 rc = reduced_cost(graph, arc, cost, potential);
-
-		assert(cap >= 0);
-		/* asserts logic implication: (rc<0 -> cap==0)*/
-		assert(!(rc < 0) || cap == 0);
-	}
-#endif
-	solved = true;
-
-finish:
-	tal_free(this_ctx);
-	return solved;
-}
-
-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)
-{
-	const tal_t *this_ctx = tal(ctx, tal_t);
-
-	assert(graph);
-	const size_t max_num_arcs = graph_max_num_arcs(graph);
-	const size_t max_num_nodes = graph_max_num_nodes(graph);
-
-	/* check preconditions */
-	assert(amount > 0);
-	assert(source.idx < max_num_nodes);
-	assert(destination.idx < max_num_nodes);
-	assert(capacity);
-	assert(cost);
-	assert(tal_count(capacity) == max_num_arcs);
-	assert(tal_count(cost) == max_num_arcs);
-
-	s64 *potential = tal_arrz(this_ctx, s64, max_num_nodes);
-	s64 *excess = tal_arrz(this_ctx, s64, max_num_nodes);
-
-	excess[source.idx] = amount;
-	excess[destination.idx] = -amount;
-
-	if (!mcf_refinement(this_ctx, graph, excess, capacity, cost, potential))
-		goto fail;
-
-	tal_free(this_ctx);
-	return true;
-
-fail:
-	tal_free(this_ctx);
-	return false;
-}
-
-s64 flow_cost(const struct graph *graph, const s64 *capacity, const s64 *cost)
-{
-	assert(graph);
-	const size_t max_num_arcs = graph_max_num_arcs(graph);
-	s64 total_cost = 0;
-
-	/* check preconditions */
-	assert(capacity);
-	assert(cost);
-	assert(tal_count(capacity) == max_num_arcs);
-	assert(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;
-}
diff --git a/plugins/askrene/algorithm.h b/plugins/askrene/algorithm.h
deleted file mode 100644
index 40010eb5..00000000
--- a/plugins/askrene/algorithm.h
+++ /dev/null
@@ -1,179 +0,0 @@
-#ifndef LIGHTNING_PLUGINS_ASKRENE_ALGORITHM_H
-#define LIGHTNING_PLUGINS_ASKRENE_ALGORITHM_H
-
-/* Implementation of network algorithms: shortests path, minimum cost flow, etc.
- */
-
-#include "config.h"
-#include <plugins/askrene/graph.h>
-
-/* Search any path from source to destination using Breadth First Search.
- *
- * input:
- * @ctx: tal allocator,
- * @graph: graph of the network,
- * @source: source node,
- * @destination: destination node,
- * @capacity: arcs capacity
- * @cap_threshold: an arc i is traversable if capacity[i]>=cap_threshold
- *
- * output:
- * @prev: prev[i] is the arc that leads to node i for an optimal solution, it
- * @return: true if the destination node was reached.
- *
- * precondition:
- * |capacity|=graph_max_num_arcs
- * |prev|=graph_max_num_nodes
- *
- * The destination is only used as a stopping condition, if destination is
- * passed with an invalid idx then the algorithm will produce a discovery tree
- * of all reacheable nodes from the source.
- * */
-bool BFS_path(const tal_t *ctx, const struct graph *graph,
-	      const struct node source, const struct node destination,
-	      const s64 *capacity, const s64 cap_threshold, struct arc *prev);
-
-
-/* Computes the distance from the source to every other node in the network
- * using Dijkstra's algorithm.
- *
- * input:
- * @ctx: tal context for internal allocation
- * @graph: topological information of the graph
- * @source: source node
- * @destination: destination node
- * @prune: if prune is true the algorithm stops when the optimal path is found
- * for the destination node
- * @capacity: arcs capacity
- * @cap_threshold: an arc i is traversable if capacity[i]>=cap_threshold
- * @cost: arc's cost
- * @potential: nodes' potential, ie. reduced cost for an arc
- * 	c_ij = cost_ij - potential[i] + potential[j]
- *
- * output:
- * @prev: for each node, this is the arc that was used to arrive to it, this can
- * be used to reconstruct the path from the destination to the source,
- * @distance: node's best distance
- * returns true if an optimal path is found for the destination, false otherwise
- *
- * precondition:
- * |capacity|=|cost|=graph_max_num_arcs
- * |prev|=|distance|=graph_max_num_nodes
- * cost[i]>=0
- * if prune is true the destination must be valid
- * */
-bool dijkstra_path(const tal_t *ctx, const struct graph *graph,
-		   const struct node source, const struct node destination,
-		   bool prune, const s64 *capacity, const s64 cap_threshold,
-		   const s64 *cost, const s64 *potential, struct arc *prev,
-		   s64 *distance);
-
-
-/* Finds any 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 simple augmenting paths algorithm.
- *
- * input:
- * @ctx: tal context for internal allocation
- * @graph: topological information of the graph
- * @source: source node
- * @destination: destination node
- * @capacity: arcs capacity
- * @amount: supply/demand
- *
- * output:
- * @capacity: residual capacity
- * returns true if the balance constraint can be satisfied
- *
- * precondition:
- * |capacity|=graph_max_num_arcs
- * amount>=0
- * */
-bool simple_feasibleflow(const tal_t *ctx, const struct graph *graph,
-			 const struct node source,
-			 const struct node destination, s64 *capacity,
-			 s64 amount);
-
-
-/* Computes the balance of a node, ie. the incoming flows minus the outgoing.
- *
- * @graph: topology
- * @node: node
- * @capacity: capacity in the residual sense, not the constrain capacity
- *
- * This works because in the adjacency list an arc wich is dual is associated
- * with an inconming arc i, then we add this flow, while an arc which is not
- * dual corresponds to and outgoing flow that we need to substract.
- * The flow on the arc i (not dual) is computed as:
- * 	flow[i] = residual_capacity[i_dual],
- * while the constrain capacity is
- * 	capacity[i] = residual_capacity[i] + residual_capacity[i_dual] */
-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);
-
-/* Take an existent flow and find an optimal redistribution:
- *
- * inputs:
- * @ctx: tal context for internal allocation,
- * @graph: topological information of the graph,
- * @excess: supply/demand of nodes,
- * @capacity: residual capacity in the arcs,
- * @cost: cost per unit of flow for every arc,
- * @potential: node potential,
- *
- * outputs:
- * @excess: all values become zero if there exist a feasible solution,
- * @capacity: encodes the resulting flow,
- * @potential: the potential that proves the solution using the complementary
- * slackness optimality condition.
- * */
-bool mcf_refinement(const tal_t *ctx,
-		    const struct graph *graph,
-		    s64 *excess,
-		    s64 *capacity,
-		    const s64 *cost,
-		    s64 *potential);
-
-#endif /* LIGHTNING_PLUGINS_ASKRENE_ALGORITHM_H */
diff --git a/plugins/askrene/askrene.c b/plugins/askrene/askrene.c
index de38b145..ff94b21a 100644
--- a/plugins/askrene/askrene.c
+++ b/plugins/askrene/askrene.c
@@ -26,9 +26,9 @@
 #include <math.h>
 #include <plugins/askrene/askrene.h>
 #include <plugins/askrene/child/child_log.h>
-#include <plugins/askrene/flow.h>
+#include <plugins/askrene/child/flow.h>
+#include <plugins/askrene/child/mcf.h>
 #include <plugins/askrene/layer.h>
-#include <plugins/askrene/mcf.h>
 #include <plugins/askrene/reserve.h>
 #include <sys/wait.h>
 #include <unistd.h>
diff --git a/plugins/askrene/child/algorithm.c b/plugins/askrene/child/algorithm.c
new file mode 100644
index 00000000..4da79e0c
--- /dev/null
+++ b/plugins/askrene/child/algorithm.c
@@ -0,0 +1,670 @@
+#include "config.h"
+#include <ccan/bitmap/bitmap.h>
+#include <ccan/tal/tal.h>
+#include <plugins/askrene/child/algorithm.h>
+#include <plugins/askrene/child/priorityqueue.h>
+
+static const s64 INFINITE = INT64_MAX;
+
+#define MAX(x, y) (((x) > (y)) ? (x) : (y))
+#define MIN(x, y) (((x) < (y)) ? (x) : (y))
+
+bool BFS_path(const tal_t *ctx, const struct graph *graph,
+	      const struct node source, const struct node destination,
+	      const s64 *capacity, const s64 cap_threshold, struct arc *prev)
+{
+	const tal_t *this_ctx = tal(ctx, tal_t);
+	bool target_found = false;
+	assert(graph);
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+	const size_t max_num_nodes = graph_max_num_nodes(graph);
+
+	/* check preconditions */
+	assert(source.idx < max_num_nodes);
+	assert(capacity);
+	assert(prev);
+	assert(tal_count(capacity) == max_num_arcs);
+	assert(tal_count(prev) == max_num_nodes);
+
+	for (size_t i = 0; i < max_num_nodes; i++)
+		prev[i].idx = INVALID_INDEX;
+
+	/* A minimalistic queue is implemented here. Nodes are not visited more
+	 * than once, therefore a maximum size of max_num_nodes is sufficient.
+	 * max_num_arcs would work as well but we expect max_num_arcs to be a
+	 * factor >10 greater than max_num_nodes. */
+	u32 *queue = tal_arr(this_ctx, u32, max_num_nodes);
+	size_t queue_start = 0, queue_end = 0;
+
+	queue[queue_end++] = source.idx;
+
+	while (queue_start < queue_end) {
+		struct node cur = {.idx = queue[queue_start++]};
+
+		if (cur.idx == destination.idx) {
+			target_found = true;
+			break;
+		}
+
+		for (struct arc arc = node_adjacency_begin(graph, cur);
+		     !node_adjacency_end(arc);
+		     arc = node_adjacency_next(graph, arc)) {
+			/* check if this arc is traversable */
+			if (capacity[arc.idx] < cap_threshold)
+				continue;
+
+			const struct node next = arc_head(graph, arc);
+
+			/* if that node has been seen previously */
+			if (prev[next.idx].idx != INVALID_INDEX ||
+			    next.idx == source.idx)
+				continue;
+
+			prev[next.idx] = arc;
+
+			assert(queue_end < max_num_nodes);
+			queue[queue_end++] = next.idx;
+		}
+	}
+
+	tal_free(this_ctx);
+	return target_found;
+}
+
+bool dijkstra_path(const tal_t *ctx, const struct graph *graph,
+		   const struct node source, const struct node destination,
+		   bool prune, const s64 *capacity, const s64 cap_threshold,
+		   const s64 *cost, const s64 *potential, struct arc *prev,
+		   s64 *distance)
+{
+	bool target_found = false;
+	assert(graph);
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+	const size_t max_num_nodes = graph_max_num_nodes(graph);
+	const tal_t *this_ctx = tal(ctx, tal_t);
+
+	/* check preconditions */
+	assert(source.idx<max_num_nodes);
+	assert(cost);
+	assert(capacity);
+	assert(prev);
+	assert(distance);
+
+	/* if prune is true then the destination cannot be invalid */
+	assert(destination.idx < max_num_nodes || !prune);
+
+	assert(tal_count(cost) == max_num_arcs);
+	assert(tal_count(capacity) == max_num_arcs);
+	assert(tal_count(prev) == max_num_nodes);
+	assert(tal_count(distance) == max_num_nodes);
+
+	/* FIXME: maybe this is unnecessary */
+	bitmap *visited = tal_arrz(this_ctx, bitmap,
+				   BITMAP_NWORDS(max_num_nodes));
+
+	for (size_t i = 0; i < max_num_nodes; ++i)
+		prev[i].idx = INVALID_INDEX;
+
+	struct priorityqueue *q;
+	q = priorityqueue_new(this_ctx, max_num_nodes);
+	const s64 *const dijkstra_distance = priorityqueue_value(q);
+
+	priorityqueue_init(q);
+	priorityqueue_update(q, source.idx, 0);
+
+	while (!priorityqueue_empty(q)) {
+		const u32 cur = priorityqueue_top(q);
+		priorityqueue_pop(q);
+
+		/* FIXME: maybe this is unnecessary */
+		if (bitmap_test_bit(visited, cur))
+			continue;
+		bitmap_set_bit(visited, cur);
+
+		if (cur == destination.idx) {
+			target_found = true;
+			if (prune)
+				break;
+		}
+
+		for (struct arc arc =
+			 node_adjacency_begin(graph, node_obj(cur));
+		     !node_adjacency_end(arc);
+		     arc = node_adjacency_next(graph, arc)) {
+			/* check if this arc is traversable */
+			if (capacity[arc.idx] < cap_threshold)
+				continue;
+
+			const struct node next = arc_head(graph, arc);
+
+			const s64 cij = cost[arc.idx] - potential[cur] +
+					potential[next.idx];
+
+			/* Dijkstra only works with non-negative weights */
+			assert(cij >= 0);
+
+			if (dijkstra_distance[next.idx] <=
+			    dijkstra_distance[cur] + cij)
+				continue;
+
+			priorityqueue_update(q, next.idx,
+					     dijkstra_distance[cur] + cij);
+			prev[next.idx] = arc;
+		}
+	}
+	for (size_t i = 0; i < max_num_nodes; i++)
+		distance[i] = dijkstra_distance[i];
+
+	tal_free(this_ctx);
+	return target_found;
+}
+
+/* Get the max amount of flow one can send from source to target along the path
+ * encoded in `prev`. */
+static s64 get_augmenting_flow(const struct graph *graph,
+			       const struct node source,
+			       const struct node target, const s64 *capacity,
+			       const struct arc *prev)
+{
+	const size_t max_num_nodes = graph_max_num_nodes(graph);
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+	assert(max_num_nodes == tal_count(prev));
+	assert(max_num_arcs == tal_count(capacity));
+
+	/* count the number of arcs in the path */
+	int path_length = 0;
+	s64 flow = INFINITE;
+
+	struct node cur = target;
+	while (cur.idx != source.idx) {
+		assert(cur.idx < max_num_nodes);
+		const struct arc arc = prev[cur.idx];
+		assert(arc.idx < max_num_arcs);
+		flow = MIN(flow, capacity[arc.idx]);
+
+		/* we are traversing in the opposite direction to the flow,
+		 * hence the next node is at the tail of the arc. */
+		cur = arc_tail(graph, arc);
+
+		/* We may never have a path exceeds the number of nodes, it this
+		 * happens it means we have an infinite loop. */
+		path_length++;
+		if(path_length >= max_num_nodes){
+			flow = -1;
+			break;
+		}
+	}
+
+	assert(flow < INFINITE && flow > 0);
+	return flow;
+}
+
+
+/* Helper.
+ * Sends an amount of flow through an arc, changing the flow balance of the
+ * nodes connected by the arc and the [residual] capacity of the arc and its
+ * dual. */
+static void sendflow(const struct graph *graph, const struct arc arc,
+		     const s64 flow, s64 *arc_capacity, s64 *node_balance)
+{
+	const struct arc dual = arc_dual(graph, arc);
+
+	arc_capacity[arc.idx] -= flow;
+	arc_capacity[dual.idx] += flow;
+
+	if (node_balance) {
+		const struct node src = arc_tail(graph, arc),
+				  dst = arc_tail(graph, dual);
+
+		node_balance[src.idx] -= flow;
+		node_balance[dst.idx] += flow;
+	}
+}
+
+/* Augment a `flow` amount along the path defined by `prev`.*/
+static void augment_flow(const struct graph *graph,
+			 const struct node source,
+			 const struct node target,
+			 const struct arc *prev,
+			 s64 *excess,
+			 s64 *capacity,
+			 s64 flow)
+{
+	const size_t max_num_nodes = graph_max_num_nodes(graph);
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+	assert(max_num_nodes == tal_count(prev));
+	assert(max_num_arcs == tal_count(capacity));
+
+	struct node cur = target;
+	/* count the number of arcs in the path */
+	int path_length = 0;
+
+	while (cur.idx != source.idx) {
+		assert(cur.idx < max_num_nodes);
+		const struct arc arc = prev[cur.idx];
+
+		sendflow(graph, arc, flow, capacity, excess);
+
+		/* we are traversing in the opposite direction to the flow,
+		 * hence the next node is at the tail of the arc. */
+		cur = arc_tail(graph, arc);
+
+		/* We may never have a path exceeds the number of nodes, it this
+		 * happens it means we have an infinite loop. */
+		path_length++;
+		if (path_length >= max_num_nodes)
+			break;
+	}
+	assert(path_length < max_num_nodes);
+}
+
+bool simple_feasibleflow(const tal_t *ctx,
+			 const struct graph *graph,
+			 const struct node source,
+			 const struct node destination,
+			 s64 *capacity,
+			 s64 amount)
+{
+	const tal_t *this_ctx = tal(ctx, tal_t);
+	assert(graph);
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+	const size_t max_num_nodes = graph_max_num_nodes(graph);
+
+	/* check preconditions */
+	assert(amount > 0);
+	assert(source.idx < max_num_nodes);
+	assert(destination.idx < max_num_nodes);
+	assert(capacity);
+	assert(tal_count(capacity) == max_num_arcs);
+
+	/* path information
+	 * prev: is the id of the arc that lead to the node. */
+	struct arc *prev = tal_arr(this_ctx, struct arc, max_num_nodes);
+	if (!prev)
+		goto finish;
+
+	while (amount > 0) {
+		/* find a path from source to target */
+		if (!BFS_path(this_ctx, graph, source, destination, capacity, 1,
+			      prev))
+			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(amount, delta);
+		assert(delta > 0 && delta <= amount);
+
+		augment_flow(graph, source, destination, prev, NULL, capacity,
+			     delta);
+		amount -= delta;
+	}
+finish:
+	tal_free(this_ctx);
+	return amount == 0;
+}
+
+s64 node_balance(const struct graph *graph,
+		 const struct node node,
+		 const s64 *capacity)
+{
+	s64 balance = 0;
+
+	for (struct arc arc = node_adjacency_begin(graph, node);
+	     !node_adjacency_end(arc); arc = node_adjacency_next(graph, arc)) {
+		struct arc dual = arc_dual(graph, arc);
+
+		if (arc_is_dual(graph, arc))
+			balance += capacity[arc.idx];
+		else
+			balance -= capacity[dual.idx];
+	}
+	return balance;
+}
+
+/* Helper.
+ * Compute the reduced cost of an arc. */
+static s64 reduced_cost(const struct graph *graph, const struct arc arc,
+			const s64 *cost, const s64 *potential)
+{
+	struct node src = arc_tail(graph, arc);
+	struct node dst = arc_head(graph, arc);
+	return cost[arc.idx] - potential[src.idx] + potential[dst.idx];
+}
+
+/* Finds an optimal path from the source to the nearest sink node, by definition
+ * a node i is a sink if node_balance[i]<0. It uses a reduced cost:
+ *	reduced_cost[i,j] = cost[i,j] - potential[i] + potential[j]
+ *
+ * */
+static struct node dijkstra_nearest_sink(const tal_t *ctx,
+					 const struct graph *graph,
+					 const struct node source,
+					 const s64 *node_balance,
+					 const s64 *capacity,
+					 const s64 cap_threshold,
+					 const s64 *cost,
+					 const s64 *potential,
+					 struct arc *prev,
+					 s64 *distance)
+{
+	struct node target = {.idx = INVALID_INDEX};
+	const tal_t *this_ctx = tal(ctx, tal_t);
+
+	/* check preconditions */
+	assert(graph);
+	assert(node_balance);
+	assert(capacity);
+	assert(cost);
+	assert(potential);
+	assert(prev);
+	assert(distance);
+
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+	const size_t max_num_nodes = graph_max_num_nodes(graph);
+
+	assert(source.idx < max_num_nodes);
+	assert(tal_count(node_balance) == max_num_nodes);
+	assert(tal_count(capacity) == max_num_arcs);
+	assert(tal_count(cost) == max_num_arcs);
+	assert(tal_count(potential) == max_num_nodes);
+	assert(tal_count(prev) == max_num_nodes);
+	assert(tal_count(distance) == max_num_nodes);
+
+	for (size_t i = 0; i < max_num_arcs; i++) {
+		/* is this arc saturated? */
+		if (capacity[i] < cap_threshold)
+			continue;
+
+		struct arc arc = {.idx = i};
+		struct node tail = arc_tail(graph, arc);
+		struct node head = arc_head(graph, arc);
+		s64 red_cost =
+		    cost[i] - potential[tail.idx] + potential[head.idx];
+
+		/* reducted cost cannot be negative for non saturated arcs,
+		 * otherwise Dijkstra does not work. */
+		if (red_cost < 0)
+			goto finish;
+	}
+
+	for (size_t i = 0; i < max_num_nodes; ++i)
+		prev[i].idx = INVALID_INDEX;
+
+/* Only in debug mode we keep track of visited nodes. */
+#ifdef ASKRENE_UNITTEST
+	bitmap *visited =
+	    tal_arrz(this_ctx, bitmap, BITMAP_NWORDS(max_num_nodes));
+#endif
+
+	struct priorityqueue *q;
+	q = priorityqueue_new(this_ctx, max_num_nodes);
+	const s64 *const dijkstra_distance = priorityqueue_value(q);
+
+	priorityqueue_init(q);
+	priorityqueue_update(q, source.idx, 0);
+
+	while (!priorityqueue_empty(q)) {
+		const u32 idx = priorityqueue_top(q);
+		const struct node cur = {.idx = idx};
+		priorityqueue_pop(q);
+
+/* Only in debug mode we keep track of visited nodes. */
+#ifdef ASKRENE_UNITTEST
+		assert(!bitmap_test_bit(visited, cur.idx));
+		bitmap_set_bit(visited, cur.idx);
+#endif
+
+		if (node_balance[cur.idx] < 0) {
+			target = cur;
+			break;
+		}
+
+		for (struct arc arc = node_adjacency_begin(graph, cur);
+		     !node_adjacency_end(arc);
+		     arc = node_adjacency_next(graph, arc)) {
+			/* check if this arc is traversable */
+			if (capacity[arc.idx] < cap_threshold)
+				continue;
+
+			const struct node next = arc_head(graph, arc);
+
+			const s64 cij = cost[arc.idx] - potential[cur.idx] +
+					potential[next.idx];
+
+			/* Dijkstra only works with non-negative weights */
+			assert(cij >= 0);
+
+			if (dijkstra_distance[next.idx] <=
+			    dijkstra_distance[cur.idx] + cij)
+				continue;
+
+			priorityqueue_update(q, next.idx,
+					     dijkstra_distance[cur.idx] + cij);
+			prev[next.idx] = arc;
+		}
+	}
+	for (size_t i = 0; i < max_num_nodes; i++)
+		distance[i] = dijkstra_distance[i];
+
+finish:
+	tal_free(this_ctx);
+	return target;
+}
+
+/* Problem: find a potential and capacity redistribution such that:
+ *	excess[all nodes] = 0
+ *	capacity[all arcs] >= 0
+ *	cost/potential [i,j] < 0 implies capacity[i,j] = 0
+ *
+ *	Q. Is this a feasible solution?
+ *
+ *	A. If we use flow conserving function sendflow, then
+ *	if for all nodes excess[i] = 0 and capacity[i,j] >= 0 for all arcs
+ *	then we have reached a feasible flow.
+ *
+ *	Q. Is this flow optimal?
+ *
+ *	A. According to Theorem 9.4 (Ahuja page 309) we have reached an optimal
+ *	solution if we are able to find a potential and flow that satisfy the
+ *	slackness optimality conditions:
+ *
+ *		if cost_reduced[i,j] > 0 then x[i,j] = 0
+ *		if 0 < x[i,j] < u[i,j] then cost_reduced[i,j] = 0
+ *		if cost_reduced[i,j] < 0 then x[i,j] = u[i,j]
+ *
+ *	In our representation the slackness optimality conditions are equivalent
+ *	to the following condition in the residual network:
+ *
+ *		cost_reduced[i,j] < 0 then capacity[i,j] = 0
+ *
+ *	Therefore yes, the solution is optimal.
+ *
+ *	Q. Why is this useful?
+ *
+ *	A. It can be used to compute a MCF from scratch or build an optimal
+ *	solution starting from a non-optimal one, eg. if we first test the
+ *	solution feasibility we already have a solution canditate, we use that
+ *	flow as input to this function, in another example we might have an
+ *	algorithm that changes the cost function at every iteration and we need
+ *	to find the MCF every time.
+ * */
+bool mcf_refinement(const tal_t *ctx,
+		    const struct graph *graph,
+		    s64 *excess,
+		    s64 *capacity,
+		    const s64 *cost,
+		    s64 *potential)
+{
+	bool solved = false;
+	const tal_t *this_ctx = tal(ctx, tal_t);
+
+	assert(graph);
+	assert(excess);
+	assert(capacity);
+	assert(cost);
+	assert(potential);
+
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+	const size_t max_num_nodes = graph_max_num_nodes(graph);
+
+	assert(tal_count(excess) == max_num_nodes);
+	assert(tal_count(capacity) == max_num_arcs);
+	assert(tal_count(cost) == max_num_arcs);
+	assert(tal_count(potential) == max_num_nodes);
+
+	s64 total_excess = 0;
+	for (u32 i = 0; i < max_num_nodes; i++)
+		total_excess += excess[i];
+
+	if (total_excess)
+		/* there is no way to satisfy the constraints if supply does not
+		 * match demand */
+		goto finish;
+
+	/* Enforce the complementary slackness condition, rolls back
+	 * constraints.  */
+	for (u32 arc_id = 0; arc_id < max_num_arcs; arc_id++) {
+		struct arc arc = {.idx = arc_id};
+		if(!arc_enabled(graph, arc))
+			continue;
+		const s64 r = capacity[arc.idx];
+		if (reduced_cost(graph, arc, cost, potential) < 0 && r > 0) {
+			/* This arc's reduced cost is negative and non
+			 * saturated. */
+			sendflow(graph, arc, r, capacity, excess);
+		}
+	}
+
+	struct arc *prev = tal_arr(this_ctx, struct arc, max_num_nodes);
+	s64 *distance = tal_arrz(this_ctx, s64, max_num_nodes);
+	if (!prev || !distance)
+		goto finish;
+
+	/* Now build back constraints again keeping the complementary slackness
+	 * condition. */
+	for (u32 node_id = 0; node_id < max_num_nodes; node_id++) {
+		struct node src = {.idx = node_id};
+
+		/* is this node a source */
+		while (excess[src.idx] > 0) {
+
+			/* where is the nearest sink */
+			struct node dst = dijkstra_nearest_sink(
+			    this_ctx, graph, src, excess, capacity, 1, cost,
+			    potential, prev, distance);
+
+			if (dst.idx >= max_num_nodes)
+				/* we failed to find a reacheable sink */
+				goto finish;
+
+			/* traverse the path and see how much flow we can send
+			 */
+			s64 delta = get_augmenting_flow(graph, src, dst,
+							capacity, prev);
+
+			delta = MIN(excess[src.idx], delta);
+			delta = MIN(-excess[dst.idx], delta);
+			assert(delta > 0);
+
+			/* commit that flow to the path */
+			augment_flow(graph, src, dst, prev, excess, capacity,
+				     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[dst.idx], distance[n]);
+			}
+		}
+	}
+
+#ifdef ASKRENE_UNITTEST
+	/* verify that we have satisfied all constraints */
+	for (u32 i = 0; i < max_num_nodes; i++) {
+		assert(excess[i] == 0);
+	}
+	for (u32 i = 0; i < max_num_arcs; i++) {
+		struct arc arc = {.idx = i};
+		if(!arc_enabled(graph, arc))
+			continue;
+		const s64 cap = capacity[arc.idx];
+		const s64 rc = reduced_cost(graph, arc, cost, potential);
+
+		assert(cap >= 0);
+		/* asserts logic implication: (rc<0 -> cap==0)*/
+		assert(!(rc < 0) || cap == 0);
+	}
+#endif
+	solved = true;
+
+finish:
+	tal_free(this_ctx);
+	return solved;
+}
+
+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)
+{
+	const tal_t *this_ctx = tal(ctx, tal_t);
+
+	assert(graph);
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+	const size_t max_num_nodes = graph_max_num_nodes(graph);
+
+	/* check preconditions */
+	assert(amount > 0);
+	assert(source.idx < max_num_nodes);
+	assert(destination.idx < max_num_nodes);
+	assert(capacity);
+	assert(cost);
+	assert(tal_count(capacity) == max_num_arcs);
+	assert(tal_count(cost) == max_num_arcs);
+
+	s64 *potential = tal_arrz(this_ctx, s64, max_num_nodes);
+	s64 *excess = tal_arrz(this_ctx, s64, max_num_nodes);
+
+	excess[source.idx] = amount;
+	excess[destination.idx] = -amount;
+
+	if (!mcf_refinement(this_ctx, graph, excess, capacity, cost, potential))
+		goto fail;
+
+	tal_free(this_ctx);
+	return true;
+
+fail:
+	tal_free(this_ctx);
+	return false;
+}
+
+s64 flow_cost(const struct graph *graph, const s64 *capacity, const s64 *cost)
+{
+	assert(graph);
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+	s64 total_cost = 0;
+
+	/* check preconditions */
+	assert(capacity);
+	assert(cost);
+	assert(tal_count(capacity) == max_num_arcs);
+	assert(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;
+}
diff --git a/plugins/askrene/child/algorithm.h b/plugins/askrene/child/algorithm.h
new file mode 100644
index 00000000..67bbff93
--- /dev/null
+++ b/plugins/askrene/child/algorithm.h
@@ -0,0 +1,179 @@
+#ifndef LIGHTNING_PLUGINS_ASKRENE_CHILD_ALGORITHM_H
+#define LIGHTNING_PLUGINS_ASKRENE_CHILD_ALGORITHM_H
+
+/* Implementation of network algorithms: shortests path, minimum cost flow, etc.
+ */
+
+#include "config.h"
+#include <plugins/askrene/child/graph.h>
+
+/* Search any path from source to destination using Breadth First Search.
+ *
+ * input:
+ * @ctx: tal allocator,
+ * @graph: graph of the network,
+ * @source: source node,
+ * @destination: destination node,
+ * @capacity: arcs capacity
+ * @cap_threshold: an arc i is traversable if capacity[i]>=cap_threshold
+ *
+ * output:
+ * @prev: prev[i] is the arc that leads to node i for an optimal solution, it
+ * @return: true if the destination node was reached.
+ *
+ * precondition:
+ * |capacity|=graph_max_num_arcs
+ * |prev|=graph_max_num_nodes
+ *
+ * The destination is only used as a stopping condition, if destination is
+ * passed with an invalid idx then the algorithm will produce a discovery tree
+ * of all reacheable nodes from the source.
+ * */
+bool BFS_path(const tal_t *ctx, const struct graph *graph,
+	      const struct node source, const struct node destination,
+	      const s64 *capacity, const s64 cap_threshold, struct arc *prev);
+
+
+/* Computes the distance from the source to every other node in the network
+ * using Dijkstra's algorithm.
+ *
+ * input:
+ * @ctx: tal context for internal allocation
+ * @graph: topological information of the graph
+ * @source: source node
+ * @destination: destination node
+ * @prune: if prune is true the algorithm stops when the optimal path is found
+ * for the destination node
+ * @capacity: arcs capacity
+ * @cap_threshold: an arc i is traversable if capacity[i]>=cap_threshold
+ * @cost: arc's cost
+ * @potential: nodes' potential, ie. reduced cost for an arc
+ * 	c_ij = cost_ij - potential[i] + potential[j]
+ *
+ * output:
+ * @prev: for each node, this is the arc that was used to arrive to it, this can
+ * be used to reconstruct the path from the destination to the source,
+ * @distance: node's best distance
+ * returns true if an optimal path is found for the destination, false otherwise
+ *
+ * precondition:
+ * |capacity|=|cost|=graph_max_num_arcs
+ * |prev|=|distance|=graph_max_num_nodes
+ * cost[i]>=0
+ * if prune is true the destination must be valid
+ * */
+bool dijkstra_path(const tal_t *ctx, const struct graph *graph,
+		   const struct node source, const struct node destination,
+		   bool prune, const s64 *capacity, const s64 cap_threshold,
+		   const s64 *cost, const s64 *potential, struct arc *prev,
+		   s64 *distance);
+
+
+/* Finds any 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 simple augmenting paths algorithm.
+ *
+ * input:
+ * @ctx: tal context for internal allocation
+ * @graph: topological information of the graph
+ * @source: source node
+ * @destination: destination node
+ * @capacity: arcs capacity
+ * @amount: supply/demand
+ *
+ * output:
+ * @capacity: residual capacity
+ * returns true if the balance constraint can be satisfied
+ *
+ * precondition:
+ * |capacity|=graph_max_num_arcs
+ * amount>=0
+ * */
+bool simple_feasibleflow(const tal_t *ctx, const struct graph *graph,
+			 const struct node source,
+			 const struct node destination, s64 *capacity,
+			 s64 amount);
+
+
+/* Computes the balance of a node, ie. the incoming flows minus the outgoing.
+ *
+ * @graph: topology
+ * @node: node
+ * @capacity: capacity in the residual sense, not the constrain capacity
+ *
+ * This works because in the adjacency list an arc wich is dual is associated
+ * with an inconming arc i, then we add this flow, while an arc which is not
+ * dual corresponds to and outgoing flow that we need to substract.
+ * The flow on the arc i (not dual) is computed as:
+ * 	flow[i] = residual_capacity[i_dual],
+ * while the constrain capacity is
+ * 	capacity[i] = residual_capacity[i] + residual_capacity[i_dual] */
+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);
+
+/* Take an existent flow and find an optimal redistribution:
+ *
+ * inputs:
+ * @ctx: tal context for internal allocation,
+ * @graph: topological information of the graph,
+ * @excess: supply/demand of nodes,
+ * @capacity: residual capacity in the arcs,
+ * @cost: cost per unit of flow for every arc,
+ * @potential: node potential,
+ *
+ * outputs:
+ * @excess: all values become zero if there exist a feasible solution,
+ * @capacity: encodes the resulting flow,
+ * @potential: the potential that proves the solution using the complementary
+ * slackness optimality condition.
+ * */
+bool mcf_refinement(const tal_t *ctx,
+		    const struct graph *graph,
+		    s64 *excess,
+		    s64 *capacity,
+		    const s64 *cost,
+		    s64 *potential);
+
+#endif /* LIGHTNING_PLUGINS_ASKRENE_CHILD_ALGORITHM_H */
diff --git a/plugins/askrene/child/dijkstra.c b/plugins/askrene/child/dijkstra.c
new file mode 100644
index 00000000..4e4ca61b
--- /dev/null
+++ b/plugins/askrene/child/dijkstra.c
@@ -0,0 +1,186 @@
+#define NDEBUG 1
+#include "config.h"
+#include <plugins/askrene/child/dijkstra.h>
+
+/* In the heap we keep node idx, but in this structure we keep the distance
+ * value associated to every node, and their position in the heap as a pointer
+ * so that we can update the nodes inside the heap when the distance label is
+ * changed.
+ *
+ * Therefore this is no longer a multipurpose heap, the node_idx must be an
+ * index between 0 and less than max_num_nodes. */
+struct dijkstra {
+	//
+	s64 *distance;
+	u32 *base;
+	u32 **heapptr;
+	size_t heapsize;
+	struct gheap_ctx gheap_ctx;
+};
+
+static const s64 INFINITE = INT64_MAX;
+
+/* Required a global dijkstra for gheap. */
+static struct dijkstra *global_dijkstra;
+
+/* The heap comparer for Dijkstra search. Since the top element must be the one
+ * with the smallest distance, we use the operator >, rather than <. */
+static int dijkstra_less_comparer(
+	const void *const ctx UNUSED,
+	const void *const a,
+	const void *const b)
+{
+	return global_dijkstra->distance[*(u32*)a]
+		> global_dijkstra->distance[*(u32*)b];
+}
+
+/* The heap move operator for Dijkstra search. */
+static void dijkstra_item_mover(void *const dst, const void *const src)
+{
+	u32 src_idx = *(u32*)src;
+	*(u32*)dst = src_idx;
+
+	// we keep track of the pointer position of each element in the heap,
+	// for easy update.
+	global_dijkstra->heapptr[src_idx] = dst;
+}
+
+/* Allocation of resources for the heap. */
+struct dijkstra *dijkstra_new(const tal_t *ctx, size_t max_num_nodes)
+{
+	struct dijkstra *dijkstra = tal(ctx, struct dijkstra);
+
+	dijkstra->distance = tal_arr(dijkstra,s64,max_num_nodes);
+	dijkstra->base = tal_arr(dijkstra,u32,max_num_nodes);
+	dijkstra->heapptr = tal_arrz(dijkstra,u32*,max_num_nodes);
+
+	dijkstra->heapsize=0;
+
+	dijkstra->gheap_ctx.fanout=2;
+	dijkstra->gheap_ctx.page_chunks=1024;
+	dijkstra->gheap_ctx.item_size=sizeof(dijkstra->base[0]);
+	dijkstra->gheap_ctx.less_comparer=dijkstra_less_comparer;
+	dijkstra->gheap_ctx.less_comparer_ctx=NULL;
+	dijkstra->gheap_ctx.item_mover=dijkstra_item_mover;
+
+	return dijkstra;
+}
+
+
+void dijkstra_init(struct dijkstra *dijkstra)
+{
+	const size_t max_num_nodes = tal_count(dijkstra->distance);
+	dijkstra->heapsize=0;
+	for(size_t i=0;i<max_num_nodes;++i)
+	{
+		dijkstra->distance[i]=INFINITE;
+		dijkstra->heapptr[i] = NULL;
+	}
+}
+size_t dijkstra_size(const struct dijkstra *dijkstra)
+{
+	return dijkstra->heapsize;
+}
+
+size_t dijkstra_maxsize(const struct dijkstra *dijkstra)
+{
+	return tal_count(dijkstra->distance);
+}
+
+static void dijkstra_append(struct dijkstra *dijkstra, u32 node_idx, s64 distance)
+{
+	assert(dijkstra_size(dijkstra) < dijkstra_maxsize(dijkstra));
+	assert(node_idx < dijkstra_maxsize(dijkstra));
+
+	const size_t pos = dijkstra->heapsize;
+
+	dijkstra->base[pos]=node_idx;
+	dijkstra->distance[node_idx]=distance;
+	dijkstra->heapptr[node_idx] = &(dijkstra->base[pos]);
+	dijkstra->heapsize++;
+}
+
+void dijkstra_update(struct dijkstra *dijkstra, u32 node_idx, s64 distance)
+{
+	assert(node_idx < dijkstra_maxsize(dijkstra));
+
+	if(!dijkstra->heapptr[node_idx])
+	{
+		// not in the heap
+		dijkstra_append(dijkstra, node_idx,distance);
+		global_dijkstra = dijkstra;
+		gheap_restore_heap_after_item_increase(
+			&dijkstra->gheap_ctx,
+			dijkstra->base,
+			dijkstra->heapsize,
+			dijkstra->heapptr[node_idx]
+				- dijkstra->base);
+		global_dijkstra = NULL;
+		return;
+	}
+
+	if(dijkstra->distance[node_idx] > distance)
+	{
+		// distance decrease
+		dijkstra->distance[node_idx] = distance;
+
+		global_dijkstra = dijkstra;
+		gheap_restore_heap_after_item_increase(
+			&dijkstra->gheap_ctx,
+			dijkstra->base,
+			dijkstra->heapsize,
+			dijkstra->heapptr[node_idx]
+				- dijkstra->base);
+		global_dijkstra = NULL;
+	}else
+	{
+		// distance increase
+		dijkstra->distance[node_idx] = distance;
+
+		global_dijkstra = dijkstra;
+		gheap_restore_heap_after_item_decrease(
+			&dijkstra->gheap_ctx,
+			dijkstra->base,
+			dijkstra->heapsize,
+			dijkstra->heapptr[node_idx]
+				- dijkstra->base);
+		global_dijkstra = NULL;
+
+	}
+  	// assert(gheap_is_heap(&dijkstra->gheap_ctx,
+	//                      dijkstra->base,
+	// 		     dijkstra_size()));
+}
+
+u32 dijkstra_top(const struct dijkstra *dijkstra)
+{
+	return dijkstra->base[0];
+}
+
+bool dijkstra_empty(const struct dijkstra *dijkstra)
+{
+	return dijkstra->heapsize==0;
+}
+
+void dijkstra_pop(struct dijkstra *dijkstra)
+{
+	if(dijkstra->heapsize==0)
+		return;
+
+	const u32 top = dijkstra_top(dijkstra);
+	assert(dijkstra->heapptr[top]==dijkstra->base);
+
+	global_dijkstra = dijkstra;
+	gheap_pop_heap(
+		&dijkstra->gheap_ctx,
+		dijkstra->base,
+		dijkstra->heapsize--);
+	global_dijkstra = NULL;
+
+	dijkstra->heapptr[top]=NULL;
+}
+
+const s64* dijkstra_distance_data(const struct dijkstra *dijkstra)
+{
+	return dijkstra->distance;
+}
diff --git a/plugins/askrene/child/dijkstra.h b/plugins/askrene/child/dijkstra.h
new file mode 100644
index 00000000..60047626
--- /dev/null
+++ b/plugins/askrene/child/dijkstra.h
@@ -0,0 +1,30 @@
+#ifndef LIGHTNING_PLUGINS_ASKRENE_CHILD_DIJKSTRA_H
+#define LIGHTNING_PLUGINS_ASKRENE_CHILD_DIJKSTRA_H
+#include "config.h"
+#include <ccan/short_types/short_types.h>
+#include <ccan/tal/tal.h>
+#include <gheap.h>
+
+/* Allocation of resources for the heap. */
+struct dijkstra *dijkstra_new(const tal_t *ctx, size_t max_num_nodes);
+
+/* Initialization of the heap for a new Dijkstra search. */
+void dijkstra_init(struct dijkstra *dijkstra);
+
+/* Inserts a new element in the heap. If node_idx was already in the heap then
+ * its distance value is updated. */
+void dijkstra_update(struct dijkstra *dijkstra, u32 node_idx, s64 distance);
+
+u32 dijkstra_top(const struct dijkstra *dijkstra);
+bool dijkstra_empty(const struct dijkstra *dijkstra);
+void dijkstra_pop(struct dijkstra *dijkstra);
+
+const s64* dijkstra_distance_data(const struct dijkstra *dijkstra);
+
+/* Number of elements on the heap. */
+size_t dijkstra_size(const struct dijkstra *dijkstra);
+
+/* Maximum number of elements the heap can host */
+size_t dijkstra_maxsize(const struct dijkstra *dijkstra);
+
+#endif /* LIGHTNING_PLUGINS_ASKRENE_CHILD_DIJKSTRA_H */
diff --git a/plugins/askrene/child/explain_failure.c b/plugins/askrene/child/explain_failure.c
new file mode 100644
index 00000000..4502a32b
--- /dev/null
+++ b/plugins/askrene/child/explain_failure.c
@@ -0,0 +1,324 @@
+#include "config.h"
+#include <ccan/tal/str/str.h>
+#include <common/dijkstra.h>
+#include <common/gossmap.h>
+#include <common/route.h>
+#include <plugins/askrene/askrene.h>
+#include <plugins/askrene/child/explain_failure.h>
+#include <plugins/askrene/layer.h>
+#include <plugins/askrene/reserve.h>
+
+#define NO_USABLE_PATHS_STRING "We could not find a usable set of paths."
+
+/* Dijkstra, reduced to ignore anything but connectivity */
+static bool always_true(const struct gossmap *map,
+			const struct gossmap_chan *c,
+			int dir,
+			struct amount_msat amount,
+			void *unused)
+{
+	return true;
+}
+
+static u64 route_score_one(struct amount_msat fee UNUSED,
+			   struct amount_msat risk UNUSED,
+			   struct amount_msat total UNUSED,
+			   int dir UNUSED,
+			   const struct gossmap_chan *c UNUSED)
+{
+	return 1;
+}
+
+/* This mirrors get_constraints() */
+static const char *why_max_constrained(const tal_t *ctx,
+				       const struct route_query *rq,
+				       struct short_channel_id_dir *scidd,
+				       struct amount_msat amount)
+{
+	char *ret = NULL;
+	const char *reservations;
+	const struct layer *constrains = NULL;
+	struct amount_msat max = amount;
+
+	/* Figure out the layer that constrains us (most) */
+	for (size_t i = 0; i < tal_count(rq->layers); i++) {
+		struct amount_msat min = AMOUNT_MSAT(0), new_max = max;
+
+		layer_apply_constraints(rq->layers[i], scidd, &min, &new_max);
+		if (!amount_msat_eq(new_max, max))
+			constrains = rq->layers[i];
+		max = new_max;
+	}
+
+	if (constrains) {
+		if (!ret)
+			ret = tal_strdup(ctx, "");
+		else
+			tal_append_fmt(&ret, ", ");
+		tal_append_fmt(&ret, "layer %s says max is %s",
+			       layer_name(constrains),
+			       fmt_amount_msat(tmpctx, max));
+	}
+
+	reservations = fmt_reservations(tmpctx, rq->reserved, scidd, rq->layers);
+	if (reservations) {
+		if (!ret)
+			ret = tal_strdup(ctx, "");
+		else
+			tal_append_fmt(&ret, " and ");
+		tal_append_fmt(&ret, "already reserved %s", reservations);
+	}
+
+	/* This seems unlikely, but don't return NULL. */
+	if (!ret)
+		ret = tal_fmt(ctx, "is constrained");
+	return ret;
+}
+
+struct stat {
+	size_t num_channels;
+	struct amount_msat capacity;
+};
+
+struct node_stats {
+	struct stat total, gossip_known, enabled;
+};
+
+enum node_direction {
+	INTO_NODE,
+	OUT_OF_NODE,
+};
+
+static void add_stat(struct stat *stat,
+		     struct amount_msat amount)
+{
+	stat->num_channels++;
+	if (!amount_msat_accumulate(&stat->capacity, amount))
+		abort();
+}
+
+static void node_stats(const struct route_query *rq,
+		       const struct gossmap_node *node,
+		       enum node_direction node_direction,
+		       struct node_stats *stats)
+{
+	memset(stats, 0, sizeof(*stats));
+	for (size_t i = 0; i < node->num_chans; i++) {
+		int dir;
+		struct gossmap_chan *c;
+		struct amount_msat cap_msat;
+
+		c = gossmap_nth_chan(rq->gossmap, node, i, &dir);
+		cap_msat = gossmap_chan_get_capacity(rq->gossmap, c);
+
+		if (node_direction == INTO_NODE)
+			dir = !dir;
+
+		add_stat(&stats->total, cap_msat);
+		if (gossmap_chan_set(c, dir))
+			add_stat(&stats->gossip_known, cap_msat);
+		if (c->half[dir].enabled)
+			add_stat(&stats->enabled, cap_msat);
+	}
+}
+
+static const char *check_capacity(const tal_t *ctx,
+				  const struct route_query *rq,
+				  const struct gossmap_node *node,
+				  enum node_direction node_direction,
+				  struct amount_msat amount,
+				  const char *name)
+{
+	struct node_stats stats;
+
+	node_stats(rq, node, node_direction, &stats);
+	if (amount_msat_greater(amount, stats.total.capacity)) {
+		return rq_log(ctx, rq, LOG_DBG,
+			      NO_USABLE_PATHS_STRING
+			      "  Total %s capacity is only %s"
+			      " (in %zu channels).",
+			      name,
+			      fmt_amount_msat(tmpctx, stats.total.capacity),
+			      stats.total.num_channels);
+	}
+	if (amount_msat_greater(amount, stats.gossip_known.capacity)) {
+		return rq_log(ctx, rq, LOG_DBG,
+			      NO_USABLE_PATHS_STRING
+			      "  Missing gossip for %s: only known %zu/%zu channels, leaving capacity only %s of %s.",
+			      name,
+			      stats.gossip_known.num_channels,
+			      stats.total.num_channels,
+			      fmt_amount_msat(tmpctx, stats.gossip_known.capacity),
+			      fmt_amount_msat(tmpctx, stats.total.capacity));
+	}
+	if (amount_msat_greater(amount, stats.enabled.capacity)) {
+		return rq_log(ctx, rq, LOG_DBG,
+			      NO_USABLE_PATHS_STRING
+			      "  The %s has disabled %zu of %zu channels, leaving capacity only %s of %s.",
+			      name,
+			      stats.total.num_channels - stats.enabled.num_channels,
+			      stats.total.num_channels,
+			      fmt_amount_msat(tmpctx, stats.enabled.capacity),
+			      fmt_amount_msat(tmpctx, stats.total.capacity));
+	}
+	return NULL;
+}
+
+/* Return description of why scidd is disabled scidd */
+static const char *describe_disabled(const tal_t *ctx,
+				     const struct route_query *rq,
+				     const struct gossmap_chan *c,
+				     const struct short_channel_id_dir *scidd)
+{
+	for (int i = tal_count(rq->layers) - 1; i >= 0; i--) {
+		struct gossmap_node *dst = gossmap_nth_node(rq->gossmap, c, !scidd->dir);
+		struct node_id dstid;
+
+		gossmap_node_get_id(rq->gossmap, dst, &dstid);
+		if (layer_disables_node(rq->layers[i], &dstid))
+			return tal_fmt(ctx, "leads to node disabled by layer %s.",
+				       layer_name(rq->layers[i]));
+		else if (layer_disables_chan(rq->layers[i], scidd)) {
+			return tal_fmt(ctx, "marked disabled by layer %s.",
+				       layer_name(rq->layers[i]));
+		}
+	}
+
+	return tal_fmt(ctx, "marked disabled by gossip message.");
+}
+
+static const char *describe_capacity(const tal_t *ctx,
+				     const struct route_query *rq,
+				     const struct short_channel_id_dir *scidd,
+				     struct amount_msat amount)
+{
+	for (int i = tal_count(rq->layers) - 1; i >= 0; i--) {
+		if (layer_created(rq->layers[i], scidd->scid)) {
+			return tal_fmt(ctx, " (created by layer %s) isn't big enough to carry %s.",
+				       layer_name(rq->layers[i]),
+				       fmt_amount_msat(tmpctx, amount));
+		}
+	}
+
+	return tal_fmt(ctx, "isn't big enough to carry %s.",
+		       fmt_amount_msat(tmpctx, amount));
+}
+
+/* We failed to find a flow at all.  Why? */
+const char *explain_failure(const tal_t *ctx,
+			    const struct route_query *rq,
+			    const struct gossmap_node *srcnode,
+			    const struct gossmap_node *dstnode,
+			    struct amount_msat amount)
+{
+	const struct route_hop *hops;
+	const struct dijkstra *dij;
+	char *path;
+	const char *cap_check;
+	const char *explanation;
+	struct short_channel_id_dir scidd;
+	struct gossmap_chan *c;
+	struct amount_msat rolling_amount;
+	struct amount_msat *path_amount;
+
+	/* Do we have enough funds? */
+	cap_check = check_capacity(ctx, rq, srcnode, OUT_OF_NODE,
+				   amount, "source");
+	if (cap_check)
+		return cap_check;
+
+	/* Does destination have enough capacity? */
+	cap_check = check_capacity(ctx, rq, dstnode, INTO_NODE,
+				   amount, "destination");
+	if (cap_check)
+		return cap_check;
+
+	/* OK, fall back to telling them why didn't shortest path
+	 * work.  This covers the "but I have a direct channel!"
+	 * case. */
+	dij = dijkstra(tmpctx, rq->gossmap, dstnode, AMOUNT_MSAT(0), 0,
+		       always_true, route_score_one, NULL);
+	hops = route_from_dijkstra(tmpctx, rq->gossmap, dij, srcnode,
+				   AMOUNT_MSAT(0), 0);
+	if (!hops)
+		return rq_log(ctx, rq, LOG_INFORM,
+			      "There is no connection between source and destination at all");
+
+	/* Description of shortest path */
+	path = tal_strdup(tmpctx, "");
+	for (size_t i = 0; i < tal_count(hops); i++) {
+		tal_append_fmt(&path, "%s%s",
+			       i > 0 ? "->" : "",
+			       fmt_short_channel_id(tmpctx, hops[i].scid));
+	}
+
+	path_amount = tal_arr(tmpctx, struct amount_msat, tal_count(hops));
+	rolling_amount = amount;
+	for (size_t i = tal_count(hops) - 1; i < tal_count(hops); i--) {
+		scidd.scid = hops[i].scid;
+		scidd.dir = hops[i].direction;
+		c = gossmap_find_chan(rq->gossmap, &scidd.scid);
+
+		path_amount[i] = rolling_amount;
+		if (!amount_msat_add_fee(&rolling_amount,
+					 c->half[scidd.dir].base_fee,
+					 c->half[scidd.dir].proportional_fee)) {
+			/* Should not happen, but since the branch exists we use
+			 * it. */
+			explanation = tal_fmt(
+			    tmpctx, "produces a fee overflow for amount %s",
+			    fmt_amount_msat(tmpctx, rolling_amount));
+			return rq_log(ctx, rq, LOG_INFORM,
+				      NO_USABLE_PATHS_STRING
+				      "  The shortest path is %s, but %s %s",
+				      path,
+				      fmt_short_channel_id_dir(tmpctx, &scidd),
+				      explanation);
+		}
+	}
+
+	/* Now walk through this: is it disabled?  Insuff capacity? */
+	for (size_t i = 0; i < tal_count(hops); i++) {
+		struct amount_msat cap_msat, min, max, htlc_max, htlc_min;
+
+		scidd.scid = hops[i].scid;
+		scidd.dir = hops[i].direction;
+		c = gossmap_find_chan(rq->gossmap, &scidd.scid);
+		cap_msat = gossmap_chan_get_capacity(rq->gossmap, c);
+		get_constraints(rq, c, scidd.dir, &min, &max);
+		htlc_max = amount_msat(fp16_to_u64(c->half[scidd.dir].htlc_max));
+		htlc_min = amount_msat(fp16_to_u64(c->half[scidd.dir].htlc_min));
+
+		if (!gossmap_chan_set(c, scidd.dir))
+			explanation = "has no gossip";
+		else if (!c->half[scidd.dir].enabled)
+			explanation = describe_disabled(tmpctx, rq, c, &scidd);
+		else if (amount_msat_greater(path_amount[i], cap_msat))
+			explanation = describe_capacity(tmpctx, rq, &scidd, path_amount[i]);
+		else if (amount_msat_greater(path_amount[i], max))
+			explanation = why_max_constrained(tmpctx, rq,
+							  &scidd, path_amount[i]);
+		else if (amount_msat_greater(path_amount[i], htlc_max))
+			explanation = tal_fmt(tmpctx,
+					      "exceeds htlc_maximum_msat ~%s",
+					      fmt_amount_msat(tmpctx, htlc_max));
+		else if (amount_msat_less(path_amount[i], htlc_min))
+			explanation = tal_fmt(tmpctx,
+					      "below htlc_minumum_msat ~%s",
+					      fmt_amount_msat(tmpctx, htlc_min));
+		else
+			continue;
+
+		return rq_log(ctx, rq, LOG_INFORM,
+			      NO_USABLE_PATHS_STRING
+			      "  The shortest path is %s, but %s %s",
+			      path,
+			      fmt_short_channel_id_dir(tmpctx, &scidd),
+			      explanation);
+	}
+
+	return rq_log(ctx, rq, LOG_BROKEN,
+		       "Actually, I'm not sure why we didn't find the"
+		       " obvious route %s: perhaps this is a bug?",
+		       path);
+}
diff --git a/plugins/askrene/child/explain_failure.h b/plugins/askrene/child/explain_failure.h
new file mode 100644
index 00000000..2914b9ee
--- /dev/null
+++ b/plugins/askrene/child/explain_failure.h
@@ -0,0 +1,16 @@
+#ifndef LIGHTNING_PLUGINS_ASKRENE_CHILD_EXPLAIN_FAILURE_H
+#define LIGHTNING_PLUGINS_ASKRENE_CHILD_EXPLAIN_FAILURE_H
+#include "config.h"
+#include <common/amount.h>
+
+struct route_query;
+struct gossmap_node;
+
+/* When MCF returns nothing, try to explain why */
+const char *explain_failure(const tal_t *ctx,
+			    const struct route_query *rq,
+			    const struct gossmap_node *srcnode,
+			    const struct gossmap_node *dstnode,
+			    struct amount_msat amount);
+
+#endif /* LIGHTNING_PLUGINS_ASKRENE_CHILD_EXPLAIN_FAILURE_H */
diff --git a/plugins/askrene/child/flow.c b/plugins/askrene/child/flow.c
new file mode 100644
index 00000000..2acd952a
--- /dev/null
+++ b/plugins/askrene/child/flow.c
@@ -0,0 +1,189 @@
+#include "config.h"
+#include <assert.h>
+#include <ccan/tal/str/str.h>
+#include <ccan/tal/tal.h>
+#include <common/fp16.h>
+#include <common/overflows.h>
+#include <math.h>
+#include <plugins/askrene/askrene.h>
+#include <plugins/askrene/child/flow.h>
+#include <plugins/libplugin.h>
+#include <stdio.h>
+
+#ifndef SUPERVERBOSE
+#define SUPERVERBOSE(...)
+#else
+#define SUPERVERBOSE_ENABLED 1
+#endif
+
+/* How much do we deliver to destination using this set of routes */
+struct amount_msat flowset_delivers(struct plugin *plugin,
+				    struct flow **flows)
+{
+	struct amount_msat final = AMOUNT_MSAT(0);
+	for (size_t i = 0; i < tal_count(flows); i++) {
+		if (!amount_msat_accumulate(&final, flows[i]->delivers)) {
+			plugin_err(plugin, "Could not add flowsat %s to %s (%zu/%zu)",
+				   fmt_amount_msat(tmpctx, flows[i]->delivers),
+				   fmt_amount_msat(tmpctx, final),
+				   i, tal_count(flows));
+		}
+	}
+	return final;
+}
+
+/* Stolen whole-cloth from @Lagrang3 in renepay's flow.c.  Wrong
+ * because of htlc overhead in reservations! */
+static double edge_probability(const struct route_query *rq,
+			       const struct short_channel_id_dir *scidd,
+			       struct amount_msat sent)
+{
+	struct amount_msat numerator, denominator;
+	struct amount_msat mincap, maxcap, additional;
+	const struct gossmap_chan *c = gossmap_find_chan(rq->gossmap, &scidd->scid);
+
+	get_constraints(rq, c, scidd->dir, &mincap, &maxcap);
+
+	/* We add an extra per-htlc reservation for the *next* HTLC, so we "over-reserve"
+	 * on local channels.  Undo that! */
+	additional = get_additional_per_htlc_cost(rq, scidd);
+	if (!amount_msat_accumulate(&mincap, additional)
+	    || !amount_msat_accumulate(&maxcap, additional))
+		abort();
+
+	if (amount_msat_less_eq(sent, mincap))
+		return 1.0;
+	else if (amount_msat_greater(sent, maxcap))
+		return 0.0;
+
+	/* Linear probability: 1 - (spend - min) / (max - min) */
+
+	/* spend > mincap, from above. */
+	if (!amount_msat_sub(&numerator, sent, mincap))
+		abort();
+	/* This can only fail is maxcap was < mincap,
+	 * so we would be captured above */
+	if (!amount_msat_sub(&denominator, maxcap, mincap))
+		abort();
+	return 1.0 - amount_msat_ratio(numerator, denominator);
+}
+
+struct amount_msat flow_spend(struct plugin *plugin, const struct flow *flow)
+{
+	const size_t pathlen = tal_count(flow->path);
+	struct amount_msat spend = flow->delivers;
+
+	for (int i = (int)pathlen - 1; i >= 0; i--) {
+		const struct half_chan *h = flow_edge(flow, i);
+		if (!amount_msat_add_fee(&spend, h->base_fee,
+					 h->proportional_fee)) {
+			plugin_err(plugin, "Could not add fee %u/%u to amount %s in %i/%zu",
+				   h->base_fee, h->proportional_fee,
+				   fmt_amount_msat(tmpctx, spend),
+				   i, pathlen);
+		}
+	}
+
+	return spend;
+}
+
+struct amount_msat flow_fee(struct plugin *plugin, const struct flow *flow)
+{
+	struct amount_msat spend = flow_spend(plugin, flow);
+	struct amount_msat fee;
+	if (!amount_msat_sub(&fee, spend, flow->delivers)) {
+		plugin_err(plugin, "Could not subtract %s from %s for fee",
+				   fmt_amount_msat(tmpctx, flow->delivers),
+				   fmt_amount_msat(tmpctx, spend));
+	}
+
+	return fee;
+}
+
+struct amount_msat flowset_fee(struct plugin *plugin, struct flow **flows)
+{
+	struct amount_msat fee = AMOUNT_MSAT(0);
+	for (size_t i = 0; i < tal_count(flows); i++) {
+		struct amount_msat this_fee = flow_fee(plugin, flows[i]);
+		if (!amount_msat_accumulate(&fee, this_fee)) {
+			plugin_err(plugin, "Could not add %s to %s for flowset fee",
+				   fmt_amount_msat(tmpctx, this_fee),
+				   fmt_amount_msat(tmpctx, fee));
+		}
+	}
+	return fee;
+}
+
+/* Helper to access the half chan at flow index idx */
+const struct half_chan *flow_edge(const struct flow *flow, size_t idx)
+{
+	assert(flow);
+	assert(idx < tal_count(flow->path));
+	return &flow->path[idx]->half[flow->dirs[idx]];
+}
+
+/* Helper function to find the success_prob for a single flow
+ *
+ * IMPORTANT: flow->success_prob is misleading, because that's the prob. of
+ * success provided that there are no other flows in the current MPP flow set.
+ * */
+double flow_probability(const struct flow *flow,
+			const struct route_query *rq)
+{
+	const size_t pathlen = tal_count(flow->path);
+	struct amount_msat spend = flow->delivers;
+	double prob = 1.0;
+
+	for (int i = (int)pathlen - 1; i >= 0; i--) {
+		const struct half_chan *h = flow_edge(flow, i);
+		struct short_channel_id_dir scidd;
+		scidd.scid = gossmap_chan_scid(rq->gossmap, flow->path[i]);
+		scidd.dir = flow->dirs[i];
+
+		prob *= edge_probability(rq, &scidd, spend);
+
+		if (!amount_msat_add_fee(&spend, h->base_fee,
+					 h->proportional_fee)) {
+			plugin_err(rq->plugin, "Could not add fee %u/%u to amount %s in %i/%zu",
+				   h->base_fee, h->proportional_fee,
+				   fmt_amount_msat(tmpctx, spend),
+				   i, pathlen);
+		}
+	}
+
+	return prob;
+}
+
+u64 flow_delay(const struct flow *flow)
+{
+	u64 delay = 0;
+	for (size_t i = 0; i < tal_count(flow->path); i++)
+		delay += flow_edge(flow, i)->delay;
+	return delay;
+}
+
+u64 flows_worst_delay(struct flow **flows)
+{
+	u64 maxdelay = 0;
+	for (size_t i = 0; i < tal_count(flows); i++) {
+		u64 delay = flow_delay(flows[i]);
+		if (delay > maxdelay)
+			maxdelay = delay;
+	}
+	return maxdelay;
+}
+
+const char *fmt_flows_step_scid(const tal_t *ctx,
+				const struct route_query *rq,
+				const struct flow *flow, size_t i)
+{
+	struct short_channel_id_dir scidd;
+
+	scidd.scid = gossmap_chan_scid(rq->gossmap, flow->path[i]);
+	scidd.dir = flow->dirs[i];
+	return fmt_short_channel_id_dir(ctx, &scidd);
+}
+
+#ifndef SUPERVERBOSE_ENABLED
+#undef SUPERVERBOSE
+#endif
diff --git a/plugins/askrene/child/flow.h b/plugins/askrene/child/flow.h
new file mode 100644
index 00000000..d424676c
--- /dev/null
+++ b/plugins/askrene/child/flow.h
@@ -0,0 +1,69 @@
+#ifndef LIGHTNING_PLUGINS_ASKRENE_CHILD_FLOW_H
+#define LIGHTNING_PLUGINS_ASKRENE_CHILD_FLOW_H
+#include "config.h"
+#include <bitcoin/short_channel_id.h>
+#include <common/amount.h>
+#include <common/gossmap.h>
+
+struct plugin;
+struct route_query;
+
+/* An actual partial flow. */
+struct flow {
+	const struct gossmap_chan **path;
+	/* The directions to traverse. */
+	int *dirs;
+	/* Amount delivered */
+	struct amount_msat delivers;
+};
+
+/* Helper to access the half chan at flow index idx */
+const struct half_chan *flow_edge(const struct flow *flow, size_t idx);
+
+/* A big number, meaning "don't bother" (not infinite, since you may add) */
+#define FLOW_INF_COST 100000000.0
+
+/* Cost function to send @f msat through @c in direction @dir,
+ * given we already have a flow of prev_flow. */
+double flow_edge_cost(const struct gossmap *gossmap,
+		      const struct gossmap_chan *c, int dir,
+		      const struct amount_msat known_min,
+		      const struct amount_msat known_max,
+		      struct amount_msat prev_flow,
+		      struct amount_msat f,
+		      double mu,
+		      double basefee_penalty,
+		      double delay_riskfactor);
+
+/* What's the success probability of this flow in isolation? */
+double flow_probability(const struct flow *flow,
+			const struct route_query *rq);
+
+/* How much do we need to send to make this flow arrive. */
+struct amount_msat flow_spend(struct plugin *plugin, const struct flow *flow);
+
+/* How much do we pay in fees to make this flow arrive. */
+struct amount_msat flow_fee(struct plugin *plugin, const struct flow *flow);
+
+/* What fee to we pay for this entire flow set? */
+struct amount_msat flowset_fee(struct plugin *plugin, struct flow **flows);
+
+/* How much does this entire flowset deliver? */
+struct amount_msat flowset_delivers(struct plugin *plugin,
+				    struct flow **flows);
+
+/* How much CLTV does this flow require? */
+u64 flow_delay(const struct flow *flow);
+
+/* Max CLTV any of these flows requires */
+u64 flows_worst_delay(struct flow **flows);
+
+const char *fmt_flows_step_scid(const tal_t *ctx,
+				const struct route_query *rq,
+				const struct flow *flow, size_t i);
+
+/* When we need to debug */
+const char *fmt_flow_full(const tal_t *ctx,
+			  const struct route_query *rq,
+			  const struct flow *flow);
+#endif /* LIGHTNING_PLUGINS_ASKRENE_CHILD_FLOW_H */
diff --git a/plugins/askrene/child/graph.c b/plugins/askrene/child/graph.c
new file mode 100644
index 00000000..7a0d9c1d
--- /dev/null
+++ b/plugins/askrene/child/graph.c
@@ -0,0 +1,65 @@
+#include "config.h"
+#include <plugins/askrene/child/graph.h>
+
+/* in the background add the actual arc or dual arc */
+static void graph_push_outbound_arc(struct graph *graph, const struct arc arc,
+				    const struct node node)
+{
+	assert(arc.idx < graph_max_num_arcs(graph));
+	assert(node.idx < graph_max_num_nodes(graph));
+
+	/* arc is already added, skip */
+	if (graph->arc_tail[arc.idx].idx != INVALID_INDEX)
+		return;
+
+	graph->arc_tail[arc.idx] = node;
+
+	const struct arc first_arc = graph->node_adjacency_first[node.idx];
+	graph->node_adjacency_next[arc.idx] = first_arc;
+	graph->node_adjacency_first[node.idx] = arc;
+}
+
+bool graph_add_arc(struct graph *graph, const struct arc arc,
+		   const struct node from, const struct node to)
+{
+	assert(from.idx < graph->max_num_nodes);
+	assert(to.idx < graph->max_num_nodes);
+
+	const struct arc dual = arc_dual(graph, arc);
+
+	if (arc.idx >= graph->max_num_arcs || dual.idx >= graph->max_num_arcs)
+		return false;
+
+	graph_push_outbound_arc(graph, arc, from);
+	graph_push_outbound_arc(graph, dual, to);
+
+	return true;
+}
+
+struct graph *graph_new(const tal_t *ctx, const size_t max_num_nodes,
+			const size_t max_num_arcs, const size_t arc_dual_bit)
+{
+	struct graph *graph;
+	graph = tal(ctx, struct graph);
+
+	graph->max_num_arcs = max_num_arcs;
+	graph->max_num_nodes = max_num_nodes;
+	graph->arc_dual_bit = arc_dual_bit;
+
+	graph->arc_tail = tal_arr(graph, struct node, graph->max_num_arcs);
+	graph->node_adjacency_first =
+	    tal_arr(graph, struct arc, graph->max_num_nodes);
+	graph->node_adjacency_next =
+	    tal_arr(graph, struct arc, graph->max_num_arcs);
+
+	/* initialize with invalid indexes so that we know these slots have
+	 * never been used, eg. arc/node is newly created */
+	for (size_t i = 0; i < graph->max_num_arcs; i++)
+		graph->arc_tail[i] = node_obj(INVALID_INDEX);
+	for (size_t i = 0; i < graph->max_num_nodes; i++)
+		graph->node_adjacency_first[i] = arc_obj(INVALID_INDEX);
+	for (size_t i = 0; i < graph->max_num_nodes; i++)
+		graph->node_adjacency_next[i] = arc_obj(INVALID_INDEX);
+
+	return graph;
+}
diff --git a/plugins/askrene/child/graph.h b/plugins/askrene/child/graph.h
new file mode 100644
index 00000000..6107fd65
--- /dev/null
+++ b/plugins/askrene/child/graph.h
@@ -0,0 +1,171 @@
+#ifndef LIGHTNING_PLUGINS_ASKRENE_CHILD_GRAPH_H
+#define LIGHTNING_PLUGINS_ASKRENE_CHILD_GRAPH_H
+
+/* Defines a graph data structure. */
+
+#include "config.h"
+#include <assert.h>
+#include <ccan/short_types/short_types.h>
+#include <ccan/tal/tal.h>
+
+#define INVALID_INDEX 0xffffffff
+
+/* A directed arc in a graph.
+ * It is a simple data object for typesafey. */
+struct arc {
+	/* arc's index */
+	u32 idx;
+};
+
+/* A node in a graph.
+ * It is a simple data object for typesafety. */
+struct node {
+	/* node's index */
+	u32 idx;
+};
+
+static inline struct arc arc_obj(u32 index)
+{
+	struct arc arc = {.idx = index};
+	return arc;
+}
+static inline struct node node_obj(u32 index)
+{
+	struct node node = {.idx = index};
+	return node;
+}
+
+/* A graph's topology. */
+struct graph {
+	/* Every arc emanates from a node, the tail.
+	 * The head of the arc is the tail of the dual. */
+	struct node *arc_tail;
+
+	/* Adjacency data for nodes. Used to move in a graph in the direction of
+	 * the arcs by looping over all arcs that exit a node.
+	 *
+	 * For every directed arc there is a dual in the opposite direction,
+	 * therefore we can use the same adjacency information to traverse in
+	 * the head to tails direction as well. */
+	struct arc *node_adjacency_next;
+	struct arc *node_adjacency_first;
+
+	size_t max_num_arcs, max_num_nodes;
+
+	/* Bit that must be flipped to obtain the dual of an arc. */
+	size_t arc_dual_bit;
+};
+
+//////////////////////////////////////////////////////////////////////////////
+
+static inline size_t graph_max_num_arcs(const struct graph *graph)
+{
+	return graph->max_num_arcs;
+}
+static inline size_t graph_max_num_nodes(const struct graph *graph)
+{
+	return graph->max_num_nodes;
+}
+
+/* Give me the dual of an arc. */
+static inline struct arc arc_dual(const struct graph *graph, struct arc arc)
+{
+	arc.idx ^= (1U << graph->arc_dual_bit);
+	return arc;
+}
+
+/* Is this arc a dual? */
+static inline bool arc_is_dual(const struct graph *graph, struct arc arc)
+{
+	return (arc.idx & (1U << graph->arc_dual_bit)) != 0;
+}
+
+/* Give me the node at the tail of an arc. */
+static inline struct node arc_tail(const struct graph *graph,
+				   const struct arc arc)
+{
+	assert(arc.idx < graph_max_num_arcs(graph));
+	return graph->arc_tail[arc.idx];
+}
+
+/* Give me the node at the head of an arc. */
+static inline struct node arc_head(const struct graph *graph,
+				   const struct arc arc)
+{
+	const struct arc dual = arc_dual(graph, arc);
+	assert(dual.idx < graph_max_num_arcs(graph));
+	return graph->arc_tail[dual.idx];
+}
+
+/* We use an arc array but not all arcs in that array do exist in the graph. */
+static inline bool arc_enabled(const struct graph *graph, const struct arc arc)
+{
+	return graph->arc_tail[arc.idx].idx < graph->max_num_nodes;
+}
+
+/* Used to loop over the arcs that exit a node.
+ *
+ * for example:
+ *
+ * void show(struct graph *graph, struct node node) {
+ * 	printf("Showing node %" PRIu32 "\n", node.idx);
+ * 	for (struct arc arc = node_adjacency_begin(graph, node);
+ * 	     !node_adjacency_end(arc);
+ * 	     arc = node_adjacency_next(graph, arc)) {
+ * 		printf("arc id: %" PRIu32 ", (%" PRIu32 " -> %" PRIu32 ")\n",
+ * 		       arc.idx,
+ * 		       arc_tail(graph, arc).idx,
+ * 		       arc_head(graph, arc).idx);
+ * 	}
+ * }
+ * */
+static inline struct arc node_adjacency_begin(const struct graph *graph,
+					      const struct node node)
+{
+	assert(node.idx < graph_max_num_nodes(graph));
+	return graph->node_adjacency_first[node.idx];
+}
+static inline bool node_adjacency_end(const struct arc arc)
+{
+	return arc.idx == INVALID_INDEX;
+}
+static inline struct arc node_adjacency_next(const struct graph *graph,
+					     const struct arc arc)
+{
+	assert(arc.idx < graph_max_num_arcs(graph));
+	return graph->node_adjacency_next[arc.idx];
+}
+
+/* Used to loop over the arcs that enter a node. */
+static inline struct arc node_rev_adjacency_begin(const struct graph *graph,
+						  const struct node node)
+{
+	return arc_dual(graph, node_adjacency_begin(graph, node));
+}
+static inline bool node_rev_adjacency_end(const struct arc arc)
+{
+	return arc.idx == INVALID_INDEX;
+}
+static inline struct arc node_rev_adjacency_next(const struct graph *graph,
+						 const struct arc arc)
+{
+	return arc_dual(graph,
+			node_adjacency_next(graph, arc_dual(graph, arc)));
+}
+
+/* This call adds an arc to the graph, it adds also the dual automatically.
+ * An arc cannot be added twice, if the caller tries to do add the same arc
+ * twice the second call is ignored.
+ * The call fails if the arc or its dual do not fit into max_num_arcs. */
+bool graph_add_arc(struct graph *graph, const struct arc arc,
+		   const struct node from, const struct node to);
+
+/* Creates a graph object. Nodes and arcs are indexed from 0 to max_num_nodes-1
+ * and max_num_arcs-1 respectively. The max_num_arcs should be big enough to
+ * accomodate also the dual arcs, ie. if the maximum index for a problem arc is
+ * I then Idual = I^(1<<arc_dual_bit) must be a valid arc index
+ * Idual<max_num_arcs. */
+struct graph *graph_new(const tal_t *ctx, const size_t max_num_nodes,
+			const size_t max_num_arcs, const size_t arc_dual_bit);
+
+#endif /* LIGHTNING_PLUGINS_ASKRENE_CHILD_GRAPH_H */
diff --git a/plugins/askrene/child/mcf.c b/plugins/askrene/child/mcf.c
new file mode 100644
index 00000000..4cca00cc
--- /dev/null
+++ b/plugins/askrene/child/mcf.c
@@ -0,0 +1,1643 @@
+#include "config.h"
+#include <assert.h>
+#include <ccan/asort/asort.h>
+#include <ccan/bitmap/bitmap.h>
+#include <ccan/list/list.h>
+#include <ccan/tal/str/str.h>
+#include <ccan/tal/tal.h>
+#include <common/utils.h>
+#include <float.h>
+#include <inttypes.h>
+#include <math.h>
+#include <plugins/askrene/askrene.h>
+#include <plugins/askrene/child/algorithm.h>
+#include <plugins/askrene/child/dijkstra.h>
+#include <plugins/askrene/child/explain_failure.h>
+#include <plugins/askrene/child/flow.h>
+#include <plugins/askrene/child/graph.h>
+#include <plugins/askrene/child/mcf.h>
+#include <plugins/askrene/child/refine.h>
+#include <stdint.h>
+
+/* # Optimal payments
+ *
+ * In this module we reduce the routing optimization problem to a linear
+ * cost optimization problem and find a solution using MCF algorithms.
+ * The optimization of the routing itself doesn't need a precise numerical
+ * solution, since we can be happy near optimal results; e.g. paying 100 msat or
+ * 101 msat for fees doesn't make any difference if we wish to deliver 1M sats.
+ * On the other hand, we are now also considering Pickhard's
+ * [1] model to improve payment reliability,
+ * hence our optimization moves to a 2D space: either we like to maximize the
+ * probability of success of a payment or minimize the routing fees, or
+ * alternatively we construct a function of the two that gives a good compromise.
+ *
+ * Therefore from now own, the definition of optimal is a matter of choice.
+ * To simplify the API of this module, we think the best way to state the
+ * problem is:
+ *
+ * 	Find a routing solution that pays the least of fees while keeping
+ * 	the probability of success above a certain value `min_probability`.
+ *
+ *
+ * # Fee Cost
+ *
+ * Routing fees is non-linear function of the payment flow x, that's true even
+ * without the base fee:
+ *
+ * 	fee_msat = base_msat + floor(millionths*x_msat / 10^6)
+ *
+ * We approximate this fee into a linear function by computing a slope `c_fee` such
+ * that:
+ *
+ * 	fee_microsat = c_fee * x_sat
+ *
+ * Function `linear_fee_cost` computes `c_fee` based on the base and
+ * proportional fees of a channel.
+ * The final product if microsat because if only
+ * the proportional fee was considered we can have c_fee = millionths.
+ * Moving to costs based in msats means we have to either truncate payments
+ * below 1ksats or estimate as 0 cost for channels with less than 1000ppm.
+ *
+ * TODO(eduardo): shall we build a linear cost function in msats?
+ *
+ * # Probability cost
+ *
+ * The probability of success P of the payment is the product of the prob. of
+ * success of forwarding parts of the payment over all routing channels. This
+ * problem is separable if we log it, and since we would like to increase P,
+ * then we can seek to minimize -log(P), and that's our prob. cost function [1].
+ *
+ * 	- log P = sum_{i} - log P_i
+ *
+ * The probability of success `P_i` of sending some flow `x` on a channel with
+ * liquidity l in the range a<=l<b is
+ *
+ * 	P_{a,b}(x) = (b-x)/(b-a); for x > a
+ * 		   = 1.         ; for x <= a
+ *
+ * Notice that unlike the similar formula in [1], the one we propose does not
+ * contain the quantization shot noise for counting states. The formula remains
+ * valid independently of the liquidity units (sats or msats).
+ *
+ * The cost associated to probability P is then -k log P, where k is some
+ * constant. For k=1 we get the following table:
+ *
+ * 	prob | cost
+ * 	-----------
+ * 	0.01 | 4.6
+ * 	0.02 | 3.9
+ * 	0.05 | 3.0
+ * 	0.10 | 2.3
+ * 	0.20 | 1.6
+ * 	0.50 | 0.69
+ * 	0.80 | 0.22
+ * 	0.90 | 0.10
+ * 	0.95 | 0.05
+ * 	0.98 | 0.02
+ * 	0.99 | 0.01
+ *
+ * Clearly -log P(x) is non-linear; we try to linearize it piecewise:
+ * split the channel into 4 arcs representing 4 liquidity regions:
+ *
+ * 	arc_0 -> [0, a)
+ * 	arc_1 -> [a, a+(b-a)*f1)
+ * 	arc_2 -> [a+(b-a)*f1, a+(b-a)*f2)
+ * 	arc_3 -> [a+(b-a)*f2, a+(b-a)*f3)
+ *
+ * where f1 = 0.5, f2 = 0.8, f3 = 0.95;
+ * We fill arc_0's capacity with complete certainty P=1, then if more flow is
+ * needed we start filling the capacity in arc_1 until the total probability
+ * of success reaches P=0.5, then arc_2 until P=1-0.8=0.2, and finally arc_3 until
+ * P=1-0.95=0.05. We don't go further than 5% prob. of success per channel.
+
+ * TODO(eduardo): this channel linearization is hard coded into
+ * `CHANNEL_PIVOTS`, maybe we can parametrize this to take values from the config file.
+ *
+ * With this choice, the slope of the linear cost function becomes:
+ *
+ * 	m_0 = 0
+ * 	m_1 = 1.38 k /(b-a)
+ * 	m_2 = 3.05 k /(b-a)
+ * 	m_3 = 9.24 k /(b-a)
+ *
+ * Notice that one of the assumptions in [2] for the MCF problem is that flows
+ * and the slope of the costs functions are integer numbers. The only way we
+ * have at hand to make it so, is to choose a universal value of `k` that scales
+ * up the slopes so that floor(m_i) is not zero for every arc.
+ *
+ * # Combine fee and prob. costs
+ *
+ * We attempt to solve the original problem of finding the solution that
+ * pays the least fees while keeping the prob. of success above a certain value,
+ * by constructing a cost function which is a linear combination of fee and
+ * prob. costs.
+ * TODO(eduardo): investigate how this procedure is justified,
+ * possibly with the use of Lagrange optimization theory.
+ *
+ * At first, prob. and fee costs live in different dimensions, they cannot be
+ * summed, it's like comparing apples and oranges.
+ * However we propose to scale the prob. cost by a global factor k that
+ * translates into the monetization of prob. cost.
+ *
+ * This was chosen empirically from examination of typical network values.
+ *
+ * # References
+ *
+ * [1] Pickhardt and Richter, https://arxiv.org/abs/2107.05322
+ * [2] R.K. Ahuja, T.L. Magnanti, and J.B. Orlin. Network Flows:
+ * Theory, Algorithms, and Applications. Prentice Hall, 1993.
+ *
+ *
+ * TODO(eduardo) it would be interesting to see:
+ * how much do we pay for reliability?
+ * Cost_fee(most reliable solution) - Cost_fee(cheapest solution)
+ *
+ * TODO(eduardo): it would be interesting to see:
+ * how likely is the most reliable path with respect to the cheapest?
+ * Prob(reliable)/Prob(cheapest) = Exp(Cost_prob(cheapest)-Cost_prob(reliable))
+ *
+ * */
+
+#define PARTS_BITS 2
+#define CHANNEL_PARTS (1 << PARTS_BITS)
+
+// These are the probability intervals we use to decompose a channel into linear
+// cost function arcs.
+static const double CHANNEL_PIVOTS[]={0,0.5,0.8,0.95};
+
+static const s64 INFINITE = INT64_MAX;
+static const s64 MU_MAX = 100;
+
+/* every payment under 1000sat will be routed through a single path */
+static const struct amount_msat SINGLE_PATH_THRESHOLD = AMOUNT_MSAT(1000000);
+
+/* Let's try this encoding of arcs:
+ * Each channel `c` has two possible directions identified by a bit
+ * `half` or `!half`, and each one of them has to be
+ * decomposed into 4 liquidity parts in order to
+ * linearize the cost function, but also to solve MCF
+ * problem we need to keep track of flows in the
+ * residual network hence we need for each directed arc
+ * in the network there must be another arc in the
+ * opposite direction refered to as it's dual. In total
+ * 1+2+1 additional bits of information:
+ *
+ * 	(chan_idx)(half)(part)(dual)
+ *
+ * That means, for each channel we need to store the
+ * information of 16 arcs. If we implement a convex-cost
+ * solver then we can reduce that number to size(half)size(dual)=4.
+ *
+ * In the adjacency of a `node` we are going to store
+ * the outgoing arcs. If we ever need to loop over the
+ * incoming arcs then we will define a reverse adjacency
+ * API.
+ * Then for each outgoing channel `(c,half)` there will
+ * be 4 parts for the actual residual capacity, hence
+ * with the dual bit set to 0:
+ *
+ * 	(c,half,0,0)
+ * 	(c,half,1,0)
+ * 	(c,half,2,0)
+ * 	(c,half,3,0)
+ *
+ * and also we need to consider the dual arcs
+ * corresponding to the channel direction `(c,!half)`
+ * (the dual has reverse direction):
+ *
+ * 	(c,!half,0,1)
+ * 	(c,!half,1,1)
+ * 	(c,!half,2,1)
+ * 	(c,!half,3,1)
+ *
+ * These are the 8 outgoing arcs relative to `node` and
+ * associated with channel `c`. The incoming arcs will
+ * be:
+ *
+ * 	(c,!half,0,0)
+ * 	(c,!half,1,0)
+ * 	(c,!half,2,0)
+ * 	(c,!half,3,0)
+ *
+ * 	(c,half,0,1)
+ * 	(c,half,1,1)
+ * 	(c,half,2,1)
+ * 	(c,half,3,1)
+ *
+ * but they will be stored as outgoing arcs on the peer
+ * node `next`.
+ *
+ * I hope this will clarify my future self when I forget.
+ *
+ * */
+
+/*
+ * We want to use the whole number here for convenience, but
+ * we can't us a union, since bit order is implementation-defined and
+ * we want chanidx on the highest bits:
+ *
+ * [ 0       1 2     3            4 5 6 ... 31 ]
+ *   dual    part    chandir      chanidx
+ */
+#define ARC_DUAL_BITOFF (0)
+#define ARC_PART_BITOFF (1)
+#define ARC_CHANDIR_BITOFF (1 + PARTS_BITS)
+#define ARC_CHANIDX_BITOFF (1 + PARTS_BITS + 1)
+#define ARC_CHANIDX_BITS  (32 - ARC_CHANIDX_BITOFF)
+
+/* How many arcs can we have for a single channel?
+ * linearization parts, both directions, and dual */
+#define ARCS_PER_CHANNEL ((size_t)1 << (PARTS_BITS + 1 + 1))
+
+static inline void arc_to_parts(struct arc arc,
+				u32 *chanidx,
+				int *chandir,
+				u32 *part,
+				bool *dual)
+{
+	if (chanidx)
+		*chanidx = (arc.idx >> ARC_CHANIDX_BITOFF);
+	if (chandir)
+		*chandir = (arc.idx >> ARC_CHANDIR_BITOFF) & 1;
+	if (part)
+		*part = (arc.idx >> ARC_PART_BITOFF) & ((1 << PARTS_BITS)-1);
+	if (dual)
+		*dual = (arc.idx >> ARC_DUAL_BITOFF) & 1;
+}
+
+static inline struct arc arc_from_parts(u32 chanidx, int chandir, u32 part, bool dual)
+{
+	struct arc arc;
+
+	assert(part < CHANNEL_PARTS);
+	assert(chandir == 0 || chandir == 1);
+	assert(chanidx < (1U << ARC_CHANIDX_BITS));
+	arc.idx = ((u32)dual << ARC_DUAL_BITOFF)
+		| (part << ARC_PART_BITOFF)
+		| ((u32)chandir << ARC_CHANDIR_BITOFF)
+		| (chanidx << ARC_CHANIDX_BITOFF);
+	return arc;
+}
+
+#define MAX(x, y) (((x) > (y)) ? (x) : (y))
+#define MIN(x, y) (((x) < (y)) ? (x) : (y))
+
+struct pay_parameters {
+	const struct route_query *rq;
+	const struct gossmap_node *source;
+	const struct gossmap_node *target;
+
+	// how much we pay
+	struct amount_msat amount;
+
+	/* base unit for computation, ie. accuracy */
+	struct amount_msat accuracy;
+
+	// channel linearization parameters
+	double cap_fraction[CHANNEL_PARTS],
+	       cost_fraction[CHANNEL_PARTS];
+
+	double delay_feefactor;
+	double base_fee_penalty;
+};
+
+/* Helper function.
+ * Given an arc of the network (not residual) give me the flow. */
+static s64 get_arc_flow(
+		const s64 *arc_residual_capacity,
+		const struct graph *graph,
+		const struct arc arc)
+{
+	assert(!arc_is_dual(graph, arc));
+	struct arc dual = arc_dual(graph, arc);
+	assert(dual.idx < tal_count(arc_residual_capacity));
+	return arc_residual_capacity[dual.idx];
+}
+
+/* Set *capacity to value, up to *cap_on_capacity.  Reduce cap_on_capacity */
+static void set_capacity(s64 *capacity, u64 value, u64 *cap_on_capacity)
+{
+	*capacity = MIN(value, *cap_on_capacity);
+	*cap_on_capacity -= *capacity;
+}
+
+/* Helper to check whether a channel is available */
+static bool channel_is_available(const struct route_query *rq,
+				 const struct gossmap_chan *chan, const int dir)
+{
+	const u32 c_idx = gossmap_chan_idx(rq->gossmap, chan);
+	return gossmap_chan_set(chan, dir) && chan->half[dir].enabled &&
+	       !bitmap_test_bit(rq->disabled_chans, c_idx * 2 + dir);
+}
+
+/* FIXME: unit test this */
+/* The probability of forwarding a payment amount given a high and low liquidity
+ * bounds.
+ * @low: the liquidity is known to be greater or equal than "low"
+ * @high: the liquidity is known to be less than "high"
+ * @amount: how much is required to forward */
+static double pickhardt_richter_probability(struct amount_msat low,
+					    struct amount_msat high,
+					    struct amount_msat amount)
+{
+	struct amount_msat all_states, good_states;
+	if (amount_msat_greater_eq(amount, high))
+		return 0.0;
+	if (!amount_msat_deduct(&amount, low))
+		return 1.0;
+	if (!amount_msat_sub(&all_states, high, low))
+		abort(); // we expect high > low
+	if (!amount_msat_sub(&good_states, all_states, amount))
+		abort(); // we expect high > amount
+	return amount_msat_ratio(good_states, all_states);
+}
+
+// TODO(eduardo): unit test this
+/* Split a directed channel into parts with linear cost function. */
+static void linearize_channel(const struct pay_parameters *params,
+			      const struct gossmap_chan *c, const int dir,
+			      s64 *capacity, double *cost)
+{
+	struct amount_msat mincap, maxcap;
+
+	/* This takes into account any payments in progress. */
+	get_constraints(params->rq, c, dir, &mincap, &maxcap);
+
+	/* Assume if min > max, min is wrong */
+	if (amount_msat_greater(mincap, maxcap))
+		mincap = maxcap;
+
+	u64 a = amount_msat_ratio_floor(mincap, params->accuracy),
+	    b = 1 + amount_msat_ratio_floor(maxcap, params->accuracy);
+
+	/* An extra bound on capacity, here we use it to reduce the flow such
+	 * that it does not exceed htlcmax.
+	 * The cap con capacity is not greater than the amount of payment units
+	 * (msat/accuracy). The way a channel is decomposed into linear cost
+	 * arcs (code below) in ascending cost order ensures that the only the
+	 * necessary capacity to forward the payment is allocated in the lower
+	 * cost arcs. This may lead to some arcs in the decomposition (at the
+	 * high cost end) to have a capacity of 0, and we can prune them while
+	 * keeping the solution optimal. */
+	u64 cap_on_capacity =
+	    MIN(amount_msat_ratio_floor(gossmap_chan_htlc_max(c, dir),
+					params->accuracy),
+		amount_msat_ratio_ceil(params->amount, params->accuracy));
+
+	set_capacity(&capacity[0], a, &cap_on_capacity);
+	cost[0]=0;
+	for(size_t i=1;i<CHANNEL_PARTS;++i)
+	{
+		set_capacity(&capacity[i], params->cap_fraction[i]*(b-a), &cap_on_capacity);
+
+		cost[i] = params->cost_fraction[i] * 1000
+			  * amount_msat_ratio(params->amount, params->accuracy)
+			  / (b - a);
+	}
+}
+
+static int cmp_u64(const u64 *a, const u64 *b, void *unused)
+{
+	if (*a < *b)
+		return -1;
+	if (*a > *b)
+		return 1;
+	return 0;
+}
+
+static int cmp_double(const double *a, const double *b, void *unused)
+{
+	if (*a < *b)
+		return -1;
+	if (*a > *b)
+		return 1;
+	return 0;
+}
+
+static double get_median_ratio(const tal_t *working_ctx,
+			       const struct graph *graph,
+			       const double *arc_prob_cost,
+			       const s64 *arc_fee_cost)
+{
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+	u64 *u64_arr = tal_arr(working_ctx, u64, max_num_arcs);
+	double *double_arr = tal_arr(working_ctx, double, max_num_arcs);
+	size_t n = 0;
+
+	for (struct arc arc = {.idx=0};arc.idx < max_num_arcs; ++arc.idx) {
+		/* scan real arcs, not unused id slots or dual arcs */
+		if (arc_is_dual(graph, arc) || !arc_enabled(graph, arc))
+			continue;
+		assert(n < max_num_arcs/2);
+		u64_arr[n] = arc_fee_cost[arc.idx];
+		double_arr[n] = arc_prob_cost[arc.idx];
+		n++;
+	}
+	asort(u64_arr, n, cmp_u64, NULL);
+	asort(double_arr, n, cmp_double, NULL);
+
+	/* Empty network, or tiny probability, nobody cares */
+	if (n == 0 || double_arr[n/2] < 0.001)
+		return 1;
+
+	/* You need to scale arc_prob_cost by this to match arc_fee_cost */
+	return u64_arr[n/2] / double_arr[n/2];
+}
+
+static void combine_cost_function(const tal_t *working_ctx,
+				  const struct graph *graph,
+				  const double *arc_prob_cost,
+				  const s64 *arc_fee_cost, const s8 *biases,
+				  s64 mu, s64 *arc_cost)
+{
+	/* probabilty and fee costs are not directly comparable!
+	 * Scale by ratio of (positive) medians. */
+	const double k =
+	    get_median_ratio(working_ctx, graph, arc_prob_cost, arc_fee_cost);
+	const double ln_30 = log(30);
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+
+	for(struct arc arc = {.idx=0};arc.idx < max_num_arcs; ++arc.idx)
+	{
+		if (arc_is_dual(graph, arc) || !arc_enabled(graph, arc))
+			continue;
+
+		const double pcost = arc_prob_cost[arc.idx];
+		const s64 fcost = arc_fee_cost[arc.idx];
+		double combined;
+		u32 chanidx;
+		int chandir;
+		s32 bias;
+
+		assert(fcost != INFINITE);
+		assert(pcost != DBL_MAX);
+		combined = fcost*mu + (MU_MAX-mu)*pcost*k;
+
+		/* Bias is in human scale, where "bigger is better" */
+		arc_to_parts(arc, &chanidx, &chandir, NULL, NULL);
+		bias = biases[(chanidx << 1) | chandir];
+		if (bias != 0) {
+			/* After some trial and error, this gives a nice
+			 * dynamic range (25 seems to be "infinite" in
+			 * practice):
+			 *    e^(-bias / (100/ln(30)))
+			 */
+			double bias_factor = exp(-bias / (100 / ln_30));
+			arc_cost[arc.idx] = combined * bias_factor;
+		} else {
+			arc_cost[arc.idx] = combined;
+		}
+		/* and the respective dual */
+		struct arc dual = arc_dual(graph, arc);
+		arc_cost[dual.idx] = -combined;
+	}
+}
+
+/* Get the fee cost associated to this directed channel.
+ * Cost is expressed as PPM of the payment.
+ *
+ * Choose and integer `c_fee` to linearize the following fee function
+ *
+ *  	fee_msat = base_msat + floor(millionths*x_msat / 10^6)
+ *
+ * into
+ *
+ *  	fee = c_fee/10^6 * x
+ *
+ *  use `base_fee_penalty` to weight the base fee and `delay_feefactor` to
+ *  weight the CLTV delay.
+ *  */
+static s64 linear_fee_cost(u32 base_fee, u32 proportional_fee, u16 cltv_delta,
+			   double base_fee_penalty,
+			   double delay_feefactor)
+{
+	s64 pfee = proportional_fee,
+	    bfee = base_fee,
+	    delay = cltv_delta;
+
+	return pfee + bfee* base_fee_penalty+ delay*delay_feefactor;
+}
+
+/* This is inversely proportional to the amount we expect to send.  Let's
+ * assume we will send ~10th of the total amount per path.  But note
+ * that it converts to parts per million! */
+static double base_fee_penalty_estimate(struct amount_msat amount)
+{
+	return amount_msat_ratio(AMOUNT_MSAT(10000000), amount);
+}
+
+static void init_linear_network(const tal_t *ctx,
+				const struct pay_parameters *params,
+				struct graph **graph, double **arc_prob_cost,
+				s64 **arc_fee_cost, s64 **arc_capacity)
+{
+	const struct gossmap *gossmap = params->rq->gossmap;
+	const size_t max_num_chans = gossmap_max_chan_idx(gossmap);
+	const size_t max_num_arcs = max_num_chans * ARCS_PER_CHANNEL;
+	const size_t max_num_nodes = gossmap_max_node_idx(gossmap);
+
+	*graph = graph_new(ctx, max_num_nodes, max_num_arcs, ARC_DUAL_BITOFF);
+	*arc_prob_cost = tal_arr(ctx, double, max_num_arcs);
+	for (size_t i = 0; i < max_num_arcs; ++i)
+		(*arc_prob_cost)[i] = DBL_MAX;
+
+	*arc_fee_cost = tal_arr(ctx, s64, max_num_arcs);
+	for (size_t i = 0; i < max_num_arcs; ++i)
+		(*arc_fee_cost)[i] = INT64_MAX;
+
+	*arc_capacity = tal_arrz(ctx, s64, max_num_arcs);
+
+	for(struct gossmap_node *node = gossmap_first_node(gossmap);
+	    node;
+	    node=gossmap_next_node(gossmap,node))
+	{
+		const u32 node_id = gossmap_node_idx(gossmap,node);
+
+		for(size_t j=0;j<node->num_chans;++j)
+		{
+			int half;
+			const struct gossmap_chan *c = gossmap_nth_chan(gossmap,
+			                                                node, j, &half);
+
+			if (!channel_is_available(params->rq, c, half))
+				continue;
+
+			/* If a channel insists on more than our total, remove it */
+			if (amount_msat_less(params->amount, gossmap_chan_htlc_min(c, half)))
+				continue;
+
+			const u32 chan_id = gossmap_chan_idx(gossmap, c);
+
+			const struct gossmap_node *next = gossmap_nth_node(gossmap,
+									   c,!half);
+
+			const u32 next_id = gossmap_node_idx(gossmap,next);
+
+			if(node_id==next_id)
+				continue;
+
+			// `cost` is the word normally used to denote cost per
+			// unit of flow in the context of MCF.
+			double prob_cost[CHANNEL_PARTS];
+			s64 capacity[CHANNEL_PARTS];
+
+			// split this channel direction to obtain the arcs
+			// that are outgoing to `node`
+			linearize_channel(params, c, half, capacity, prob_cost);
+
+			/* linear fee_cost per unit of flow */
+			const s64 fee_cost = linear_fee_cost(
+				c->half[half].base_fee,
+				c->half[half].proportional_fee,
+				c->half[half].delay,
+				params->base_fee_penalty,
+				params->delay_feefactor);
+
+			// let's subscribe the 4 parts of the channel direction
+			// (c,half), the dual of these guys will be subscribed
+			// when the `i` hits the `next` node.
+			for(size_t k=0;k<CHANNEL_PARTS;++k)
+			{
+				/* prune arcs with 0 capacity */
+				if (capacity[k] == 0)
+					continue;
+
+				struct arc arc = arc_from_parts(chan_id, half, k, false);
+
+				graph_add_arc(*graph, arc,
+					      node_obj(node_id),
+					      node_obj(next_id));
+
+				(*arc_capacity)[arc.idx] = capacity[k];
+				(*arc_prob_cost)[arc.idx] = prob_cost[k];
+				(*arc_fee_cost)[arc.idx] = fee_cost;
+
+				// + the respective dual
+				struct arc dual = arc_dual(*graph, arc);
+
+				(*arc_capacity)[dual.idx] = 0;
+				(*arc_prob_cost)[dual.idx] = -prob_cost[k];
+				(*arc_fee_cost)[dual.idx] = -fee_cost;
+			}
+		}
+	}
+}
+
+// flow on directed channels
+struct chan_flow
+{
+	s64 half[2];
+};
+
+/* Search in the network a path of positive flow until we reach a node with
+ * positive balance (returns a node idx with positive balance)
+ * or we discover a cycle (returns a node idx with 0 balance).
+ * */
+static struct node find_path_or_cycle(
+		const tal_t *working_ctx,
+		const struct route_query *rq,
+		const struct chan_flow *chan_flow,
+		const struct node source,
+		const s64 *balance,
+
+		const struct gossmap_chan **prev_chan,
+		int *prev_dir,
+		u32 *prev_idx)
+{
+	const struct gossmap *gossmap = rq->gossmap;
+	const size_t max_num_nodes = gossmap_max_node_idx(gossmap);
+	bitmap *visited =
+	    tal_arrz(working_ctx, bitmap, BITMAP_NWORDS(max_num_nodes));
+	u32 final_idx = source.idx;
+	bitmap_set_bit(visited, final_idx);
+
+	/* It is guaranteed to halt, because we either find a node with
+	 * balance[]>0 or we hit a node twice and we stop. */
+	while (balance[final_idx] <= 0) {
+		u32 updated_idx = INVALID_INDEX;
+		struct gossmap_node *cur =
+		    gossmap_node_byidx(gossmap, final_idx);
+
+		for (size_t i = 0; i < cur->num_chans; ++i) {
+			int dir;
+			const struct gossmap_chan *c =
+			    gossmap_nth_chan(gossmap, cur, i, &dir);
+
+			if (!channel_is_available(rq, c, dir))
+				continue;
+
+			const u32 c_idx = gossmap_chan_idx(gossmap, c);
+
+			/* follow the flow */
+			if (chan_flow[c_idx].half[dir] > 0) {
+				const struct gossmap_node *n =
+				    gossmap_nth_node(gossmap, c, !dir);
+				u32 next_idx = gossmap_node_idx(gossmap, n);
+
+				prev_dir[next_idx] = dir;
+				prev_chan[next_idx] = c;
+				prev_idx[next_idx] = final_idx;
+
+				updated_idx = next_idx;
+				break;
+			}
+		}
+
+		assert(updated_idx != INVALID_INDEX);
+		assert(updated_idx != final_idx);
+		final_idx = updated_idx;
+
+		if (bitmap_test_bit(visited, updated_idx)) {
+			/* We have seen this node before, we've found a cycle.
+			 */
+			assert(balance[updated_idx] <= 0);
+			break;
+		}
+		bitmap_set_bit(visited, updated_idx);
+	}
+	return node_obj(final_idx);
+}
+
+struct list_data
+{
+	struct list_node list;
+	struct flow *flow_path;
+};
+
+/* Given a path from a node with negative balance to a node with positive
+ * balance, compute the bigest flow and substract it from the nodes balance and
+ * the channels allocation. */
+static struct flow *substract_flow(const tal_t *ctx,
+				   const struct pay_parameters *params,
+				   const struct node source,
+				   const struct node sink,
+				   s64 *balance, struct chan_flow *chan_flow,
+				   const u32 *prev_idx, const int *prev_dir,
+				   const struct gossmap_chan *const *prev_chan)
+{
+	const struct gossmap *gossmap = params->rq->gossmap;
+	assert(balance[source.idx] < 0);
+	assert(balance[sink.idx] > 0);
+	s64 delta = -balance[source.idx];
+	size_t length = 0;
+	delta = MIN(delta, balance[sink.idx]);
+
+	/* We can only walk backwards, now get me the legth of the path and the
+	 * max flow we can send through this route. */
+	for (u32 cur_idx = sink.idx; cur_idx != source.idx;
+	     cur_idx = prev_idx[cur_idx]) {
+		assert(cur_idx != INVALID_INDEX);
+		const int dir = prev_dir[cur_idx];
+		const struct gossmap_chan *const chan = prev_chan[cur_idx];
+
+		/* we could optimize here by caching the idx of the channels in
+		 * the path, but the bottleneck of the algorithm is the MCF
+		 * computation not here. */
+		const u32 chan_idx = gossmap_chan_idx(gossmap, chan);
+
+		delta = MIN(delta, chan_flow[chan_idx].half[dir]);
+		length++;
+	}
+
+	struct flow *f = tal(ctx, struct flow);
+	f->path = tal_arr(f, const struct gossmap_chan *, length);
+	f->dirs = tal_arr(f, int, length);
+
+	/* Walk again and substract the flow value (delta). */
+	assert(delta > 0);
+	balance[source.idx] += delta;
+	balance[sink.idx] -= delta;
+	for (u32 cur_idx = sink.idx; cur_idx != source.idx;
+	     cur_idx = prev_idx[cur_idx]) {
+		const int dir = prev_dir[cur_idx];
+		const struct gossmap_chan *const chan = prev_chan[cur_idx];
+		const u32 chan_idx = gossmap_chan_idx(gossmap, chan);
+
+		length--;
+		/* f->path and f->dirs contain the channels in the path in the
+		 * correct order. */
+		f->path[length] = chan;
+		f->dirs[length] = dir;
+
+		chan_flow[chan_idx].half[dir] -= delta;
+	}
+	if (!amount_msat_mul(&f->delivers, params->accuracy, delta))
+		abort();
+	return f;
+}
+
+/* Substract a flow cycle from the channel allocation. */
+static void substract_cycle(const struct gossmap *gossmap,
+			    const struct node sink,
+			    struct chan_flow *chan_flow, const u32 *prev_idx,
+			    const int *prev_dir,
+			    const struct gossmap_chan *const *prev_chan)
+{
+	s64 delta = INFINITE;
+	u32 cur_idx;
+
+	/* Compute greatest flow in this cycle. */
+	for (cur_idx = sink.idx; cur_idx!=INVALID_INDEX;) {
+		const int dir = prev_dir[cur_idx];
+		const struct gossmap_chan *const chan = prev_chan[cur_idx];
+		const u32 chan_idx = gossmap_chan_idx(gossmap, chan);
+
+		delta = MIN(delta, chan_flow[chan_idx].half[dir]);
+
+		cur_idx = prev_idx[cur_idx];
+		if (cur_idx == sink.idx)
+			/* we have come back full circle */
+			break;
+	}
+	assert(cur_idx==sink.idx);
+
+	/* Walk again and substract the flow value (delta). */
+	assert(delta < INFINITE);
+	assert(delta > 0);
+
+	for (cur_idx = sink.idx;cur_idx!=INVALID_INDEX;) {
+		const int dir = prev_dir[cur_idx];
+		const struct gossmap_chan *const chan = prev_chan[cur_idx];
+		const u32 chan_idx = gossmap_chan_idx(gossmap, chan);
+
+		chan_flow[chan_idx].half[dir] -= delta;
+
+		cur_idx = prev_idx[cur_idx];
+		if (cur_idx == sink.idx)
+			/* we have come back full circle */
+			break;
+	}
+	assert(cur_idx==sink.idx);
+}
+
+/* Given a flow in the residual network, build a set of payment flows in the
+ * gossmap that corresponds to this flow. */
+static struct flow **
+get_flow_paths(const tal_t *ctx,
+	       const tal_t *working_ctx,
+	       const struct pay_parameters *params,
+	       const struct graph *graph,
+	       const s64 *arc_residual_capacity)
+{
+	struct flow **flows = tal_arr(ctx,struct flow*,0);
+
+	const size_t max_num_chans = gossmap_max_chan_idx(params->rq->gossmap);
+	struct chan_flow *chan_flow = tal_arrz(working_ctx,struct chan_flow,max_num_chans);
+
+	const size_t max_num_nodes = gossmap_max_node_idx(params->rq->gossmap);
+	s64 *balance = tal_arrz(working_ctx,s64,max_num_nodes);
+
+	const struct gossmap_chan **prev_chan
+		= tal_arr(working_ctx,const struct gossmap_chan *,max_num_nodes);
+
+
+	int *prev_dir = tal_arr(working_ctx,int,max_num_nodes);
+	u32 *prev_idx = tal_arr(working_ctx, u32, max_num_nodes);
+
+	for (u32 node_idx = 0; node_idx < max_num_nodes; node_idx++)
+		prev_idx[node_idx] = INVALID_INDEX;
+
+	// Convert the arc based residual network flow into a flow in the
+	// directed channel network.
+	// Compute balance on the nodes.
+	for (struct node n = {.idx = 0}; n.idx < max_num_nodes; n.idx++) {
+		for(struct arc arc = node_adjacency_begin(graph,n);
+		        !node_adjacency_end(arc);
+			arc = node_adjacency_next(graph,arc))
+		{
+			if(arc_is_dual(graph, arc))
+				continue;
+			struct node m = arc_head(graph,arc);
+			s64 flow = get_arc_flow(arc_residual_capacity,
+						graph, arc);
+			u32 chanidx;
+			int chandir;
+
+			balance[n.idx] -= flow;
+			balance[m.idx] += flow;
+
+			arc_to_parts(arc, &chanidx, &chandir, NULL, NULL);
+			chan_flow[chanidx].half[chandir] +=flow;
+		}
+	}
+
+	// Select all nodes with negative balance and find a flow that reaches a
+	// positive balance node.
+	for (struct node source = {.idx = 0}; source.idx < max_num_nodes;
+	     source.idx++) {
+		// this node has negative balance, flows leaves from here
+		while (balance[source.idx] < 0) {
+			prev_chan[source.idx] = NULL;
+			struct node sink = find_path_or_cycle(
+			    working_ctx, params->rq, chan_flow, source,
+			    balance, prev_chan, prev_dir, prev_idx);
+
+			if (balance[sink.idx] > 0)
+			/* case 1. found a path */
+			{
+				struct flow *fp = substract_flow(
+				    flows, params, source, sink, balance,
+				    chan_flow, prev_idx, prev_dir, prev_chan);
+
+				tal_arr_expand(&flows, fp);
+			} else
+			/* case 2. found a cycle */
+			{
+				substract_cycle(params->rq->gossmap, sink, chan_flow,
+						prev_idx, prev_dir, prev_chan);
+			}
+		}
+	}
+	return flows;
+}
+
+/* Given a single path build a flow set. */
+static struct flow **
+get_flow_singlepath(const tal_t *ctx, const struct pay_parameters *params,
+		    const struct graph *graph, const struct gossmap *gossmap,
+		    const struct node source, const struct node destination,
+		    const u64 pay_amount, const struct arc *prev)
+{
+	struct flow **flows, *f;
+	flows = tal_arr(ctx, struct flow *, 1);
+	f = flows[0] = tal(flows, struct flow);
+
+	size_t length = 0;
+
+	for (u32 cur_idx = destination.idx; cur_idx != source.idx;) {
+		assert(cur_idx != INVALID_INDEX);
+		length++;
+		struct arc arc = prev[cur_idx];
+		struct node next = arc_tail(graph, arc);
+		cur_idx = next.idx;
+	}
+	f->path = tal_arr(f, const struct gossmap_chan *, length);
+	f->dirs = tal_arr(f, int, length);
+
+	for (u32 cur_idx = destination.idx; cur_idx != source.idx;) {
+		int chandir;
+		u32 chanidx;
+		struct arc arc = prev[cur_idx];
+		arc_to_parts(arc, &chanidx, &chandir, NULL, NULL);
+
+		length--;
+		f->path[length] = gossmap_chan_byidx(gossmap, chanidx);
+		f->dirs[length] = chandir;
+
+		struct node next = arc_tail(graph, arc);
+		cur_idx = next.idx;
+	}
+	f->delivers = params->amount;
+	return flows;
+}
+
+// TODO(eduardo): choose some default values for the minflow parameters
+/* eduardo: I think it should be clear that this module deals with linear
+ * flows, ie. base fees are not considered. Hence a flow along a path is
+ * described with a sequence of directed channels and one amount.
+ * In the `pay_flow` module there are dedicated routes to compute the actual
+ * amount to be forward on each hop.
+ *
+ * TODO(eduardo): notice that we don't pay fees to forward payments with local
+ * channels and we can tell with absolute certainty the liquidity on them.
+ * Check that local channels have fee costs = 0 and bounds with certainty (min=max). */
+// TODO(eduardo): we should LOG_DBG the process of finding the MCF while
+// adjusting the frugality factor.
+static struct flow **minflow(const tal_t *ctx,
+			     const struct route_query *rq,
+			     const struct gossmap_node *source,
+			     const struct gossmap_node *target,
+			     struct amount_msat amount,
+			     u32 mu,
+			     double delay_feefactor)
+{
+	struct flow **flow_paths;
+	/* We allocate everything off this, and free it at the end,
+	 * as we can be called multiple times without cleaning tmpctx! */
+	tal_t *working_ctx = tal(NULL, char);
+	struct pay_parameters *params = tal(working_ctx, struct pay_parameters);
+
+	params->rq = rq;
+	params->source = source;
+	params->target = target;
+	params->amount = amount;
+	/* -> We reduce the granularity of the flow by limiting the subdivision
+	 * of the payment amount into 1000 units of flow. That reduces the
+	 * computational burden for algorithms that depend on it, eg. "capacity
+	 * scaling" and "successive shortest path".
+	 * -> Using Ceil operation instead of Floor so that
+	 *      accuracy x 1000 >= amount
+	 * */
+	params->accuracy =
+	    amount_msat_max(AMOUNT_MSAT(1), amount_msat_div_ceil(amount, 1000));
+
+	// template the channel partition into linear arcs
+	params->cap_fraction[0]=0;
+	params->cost_fraction[0]=0;
+	for(size_t i =1;i<CHANNEL_PARTS;++i)
+	{
+		params->cap_fraction[i]=CHANNEL_PIVOTS[i]-CHANNEL_PIVOTS[i-1];
+		params->cost_fraction[i]=
+			log((1-CHANNEL_PIVOTS[i-1])/(1-CHANNEL_PIVOTS[i]))
+			/params->cap_fraction[i];
+	}
+
+	params->delay_feefactor = delay_feefactor;
+	params->base_fee_penalty = base_fee_penalty_estimate(amount);
+
+	// build the uncertainty network with linearization and residual arcs
+	struct graph *graph;
+	double *arc_prob_cost;
+	s64 *arc_fee_cost;
+	s64 *arc_capacity;
+	init_linear_network(working_ctx, params, &graph, &arc_prob_cost,
+			    &arc_fee_cost, &arc_capacity);
+
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+	const size_t max_num_nodes = graph_max_num_nodes(graph);
+	s64 *arc_cost;
+	s64 *node_potential;
+	s64 *node_excess;
+	arc_cost = tal_arrz(working_ctx, s64, max_num_arcs);
+	node_potential = tal_arrz(working_ctx, s64, max_num_nodes);
+	node_excess = tal_arrz(working_ctx, s64, max_num_nodes);
+
+	const struct node dst = {.idx = gossmap_node_idx(rq->gossmap, target)};
+	const struct node src = {.idx = gossmap_node_idx(rq->gossmap, source)};
+
+
+	/* Since we have constraint accuracy, ask to find a payment solution
+	 * that can pay a bit more than the actual value rathen than undershoot it.
+	 * That's why we use the ceil function here. */
+	const u64 pay_amount =
+	    amount_msat_ratio_ceil(params->amount, params->accuracy);
+
+	if (!simple_feasibleflow(working_ctx, graph, src, dst,
+				 arc_capacity, pay_amount)) {
+		rq_log(tmpctx, rq, LOG_INFORM,
+		       "%s failed: unable to find a feasible flow.", __func__);
+		goto fail;
+	}
+	combine_cost_function(working_ctx, graph, arc_prob_cost, arc_fee_cost,
+			      rq->biases, mu, arc_cost);
+
+	/* We solve a linear MCF problem. */
+	if (!mcf_refinement(working_ctx,
+			    graph,
+			    node_excess,
+			    arc_capacity,
+			    arc_cost,
+			    node_potential)) {
+		rq_log(tmpctx, rq, LOG_BROKEN,
+		       "%s: MCF optimization step failed", __func__);
+		goto fail;
+	}
+
+	/* We dissect the solution of the MCF into payment routes.
+	 * Actual amounts considering fees are computed for every
+	 * channel in the routes. */
+	flow_paths = get_flow_paths(ctx, working_ctx, params,
+				    graph, arc_capacity);
+	if(!flow_paths){
+		rq_log(tmpctx, rq, LOG_BROKEN,
+		       "%s: failed to extract flow paths from the MCF solution",
+		       __func__);
+		goto fail;
+	}
+	tal_free(working_ctx);
+	return flow_paths;
+
+fail:
+	tal_free(working_ctx);
+	return NULL;
+}
+
+/* Initialize the data vectors for the single-path solver. */
+static void init_linear_network_single_path(
+    const tal_t *ctx, const struct pay_parameters *params, struct graph **graph,
+    double **arc_prob_cost, s64 **arc_fee_cost, s64 **arc_capacity)
+{
+	const size_t max_num_chans = gossmap_max_chan_idx(params->rq->gossmap);
+	const size_t max_num_arcs = max_num_chans * ARCS_PER_CHANNEL;
+	const size_t max_num_nodes = gossmap_max_node_idx(params->rq->gossmap);
+
+	*graph = graph_new(ctx, max_num_nodes, max_num_arcs, ARC_DUAL_BITOFF);
+	*arc_prob_cost = tal_arr(ctx, double, max_num_arcs);
+	for (size_t i = 0; i < max_num_arcs; ++i)
+		(*arc_prob_cost)[i] = DBL_MAX;
+
+	*arc_fee_cost = tal_arr(ctx, s64, max_num_arcs);
+	for (size_t i = 0; i < max_num_arcs; ++i)
+		(*arc_fee_cost)[i] = INT64_MAX;
+	*arc_capacity = tal_arrz(ctx, s64, max_num_arcs);
+
+	const struct gossmap *gossmap = params->rq->gossmap;
+
+	for (struct gossmap_node *node = gossmap_first_node(gossmap); node;
+	     node = gossmap_next_node(gossmap, node)) {
+		const u32 node_id = gossmap_node_idx(gossmap, node);
+
+		for (size_t j = 0; j < node->num_chans; ++j) {
+			int half;
+			const struct gossmap_chan *c =
+			    gossmap_nth_chan(gossmap, node, j, &half);
+                        struct amount_msat mincap, maxcap;
+
+			if (!channel_is_available(params->rq, c, half))
+				continue;
+
+			/* If a channel cannot forward the total amount we don't
+			 * use it. */
+			if (amount_msat_less(params->amount,
+					     gossmap_chan_htlc_min(c, half)) ||
+			    amount_msat_greater(params->amount,
+						gossmap_chan_htlc_max(c, half)))
+				continue;
+
+			get_constraints(params->rq, c, half, &mincap, &maxcap);
+			/* Assume if min > max, min is wrong */
+			if (amount_msat_greater(mincap, maxcap))
+				mincap = maxcap;
+			/* It is preferable to work on 1msat past the known
+			 * bound. */
+			if (!amount_msat_accumulate(&maxcap, amount_msat(1)))
+				abort();
+
+			/* If amount is greater than the known liquidity upper
+			 * bound we get infinite probability cost. */
+			if (amount_msat_greater_eq(params->amount, maxcap))
+				continue;
+
+			const u32 chan_id = gossmap_chan_idx(gossmap, c);
+
+			const struct gossmap_node *next =
+			    gossmap_nth_node(gossmap, c, !half);
+
+			const u32 next_id = gossmap_node_idx(gossmap, next);
+
+			/* channel to self? */
+			if (node_id == next_id)
+				continue;
+
+			struct arc arc =
+			    arc_from_parts(chan_id, half, 0, false);
+
+			graph_add_arc(*graph, arc, node_obj(node_id),
+				      node_obj(next_id));
+
+			(*arc_capacity)[arc.idx] = 1;
+			(*arc_prob_cost)[arc.idx] =
+			    (-1.0) * log(pickhardt_richter_probability(
+					 mincap, maxcap, params->amount));
+
+			struct amount_msat fee;
+			if (!amount_msat_fee(&fee, params->amount,
+					     c->half[half].base_fee,
+					     c->half[half].proportional_fee))
+				abort();
+			(*arc_fee_cost)[arc.idx] =
+			    fee.millisatoshis + /* Raw: fee cost */
+			    params->delay_feefactor * c->half[half].delay;
+		}
+	}
+}
+
+/**
+ * API for min cost single path.
+ * @ctx: context to allocate returned flows from
+ * @rq: the route_query we're processing (for logging)
+ * @source: the source to start from
+ * @target: the target to pay
+ * @amount: the amount we want to reach @target
+ * @mu: 0 = corresponds to only probabilities, 100 corresponds to only fee.
+ * @delay_feefactor: convert 1 block delay into msat.
+ *
+ * @delay_feefactor converts 1 block delay into msat, as if it were an additional
+ * fee.  So if a CLTV delay on a node is 5 blocks, that's treated as if it
+ * were a fee of 5 * @delay_feefactor.
+ *
+ * Returns an array with one flow which deliver amount to target, or NULL.
+ */
+static struct flow **single_path_flow(const tal_t *ctx, const struct route_query *rq,
+				      const struct gossmap_node *source,
+				      const struct gossmap_node *target,
+				      struct amount_msat amount, u32 mu,
+				      double delay_feefactor)
+{
+	struct flow **flow_paths;
+	/* We allocate everything off this, and free it at the end,
+	 * as we can be called multiple times without cleaning tmpctx! */
+	tal_t *working_ctx = tal(NULL, char);
+	struct pay_parameters *params = tal(working_ctx, struct pay_parameters);
+
+	params->rq = rq;
+	params->source = source;
+	params->target = target;
+	params->amount = amount;
+	/* for the single-path solver the accuracy does not detriment
+	 * performance */
+	params->accuracy = amount;
+	params->delay_feefactor = delay_feefactor;
+	params->base_fee_penalty = base_fee_penalty_estimate(amount);
+
+	struct graph *graph;
+	double *arc_prob_cost;
+	s64 *arc_fee_cost;
+	s64 *arc_capacity;
+
+	init_linear_network_single_path(working_ctx, params, &graph,
+					&arc_prob_cost, &arc_fee_cost,
+					&arc_capacity);
+
+	const struct node dst = {.idx = gossmap_node_idx(rq->gossmap, target)};
+	const struct node src = {.idx = gossmap_node_idx(rq->gossmap, source)};
+
+	const size_t max_num_nodes = graph_max_num_nodes(graph);
+	const size_t max_num_arcs = graph_max_num_arcs(graph);
+
+	s64 *potential = tal_arrz(working_ctx, s64, max_num_nodes);
+	s64 *distance = tal_arrz(working_ctx, s64, max_num_nodes);
+	s64 *arc_cost = tal_arrz(working_ctx, s64, max_num_arcs);
+	struct arc *prev = tal_arrz(working_ctx, struct arc, max_num_nodes);
+
+	combine_cost_function(working_ctx, graph, arc_prob_cost, arc_fee_cost,
+			      rq->biases, mu, arc_cost);
+
+	/* We solve a linear cost flow problem. */
+	if (!dijkstra_path(working_ctx, graph, src, dst,
+			   /* prune = */ true, arc_capacity,
+			   /*threshold = */ 1, arc_cost, potential, prev,
+			   distance)) {
+                /* This might fail if we are unable to find a suitable route, it
+                 * doesn't mean the plugin is broken, that's why we LOG_INFORM. */
+		rq_log(tmpctx, rq, LOG_INFORM,
+		       "%s: could not find a feasible single path", __func__);
+		goto fail;
+	}
+	const u64 pay_amount =
+	    amount_msat_ratio_ceil(params->amount, params->accuracy);
+
+	/* We dissect the flow into payment routes.
+	 * Actual amounts considering fees are computed for every
+	 * channel in the routes. */
+	flow_paths = get_flow_singlepath(ctx, params, graph, rq->gossmap,
+					 src, dst, pay_amount, prev);
+	if (!flow_paths) {
+		rq_log(tmpctx, rq, LOG_BROKEN,
+		       "%s: failed to extract flow paths from the single-path "
+		       "solution",
+		       __func__);
+		goto fail;
+	}
+	if (tal_count(flow_paths) != 1) {
+		rq_log(
+		    tmpctx, rq, LOG_BROKEN,
+		    "%s: single-path solution returned a multi route solution",
+		    __func__);
+		goto fail;
+	}
+	tal_free(working_ctx);
+	return flow_paths;
+
+fail:
+	tal_free(working_ctx);
+	return NULL;
+}
+
+/* Get the scidd for the i'th hop in flow */
+static void get_scidd(const struct gossmap *gossmap, const struct flow *flow,
+		      size_t i, struct short_channel_id_dir *scidd)
+{
+	scidd->scid = gossmap_chan_scid(gossmap, flow->path[i]);
+	scidd->dir = flow->dirs[i];
+}
+
+/* We use an fp16_t approximatin for htlc_max/min: this gets the exact value. */
+static struct amount_msat
+get_chan_htlc_max(const struct route_query *rq, const struct gossmap_chan *c,
+		  const struct short_channel_id_dir *scidd)
+{
+	struct amount_msat htlc_max;
+
+	gossmap_chan_get_update_details(rq->gossmap, c, scidd->dir, NULL, NULL,
+					NULL, NULL, NULL, NULL, NULL,
+					&htlc_max);
+	return htlc_max;
+}
+
+static struct amount_msat
+get_chan_htlc_min(const struct route_query *rq, const struct gossmap_chan *c,
+		  const struct short_channel_id_dir *scidd)
+{
+	struct amount_msat htlc_min;
+
+	gossmap_chan_get_update_details(rq->gossmap, c, scidd->dir, NULL, NULL,
+					NULL, NULL, NULL, NULL, &htlc_min,
+					NULL);
+	return htlc_min;
+}
+
+static bool check_htlc_min_limits(struct route_query *rq, struct flow **flows)
+{
+
+	for (size_t k = 0; k < tal_count(flows); k++) {
+		struct flow *flow = flows[k];
+		size_t pathlen = tal_count(flow->path);
+		struct amount_msat hop_amt = flow->delivers;
+		for (size_t i = pathlen - 1; i < pathlen; i--) {
+			const struct half_chan *h = flow_edge(flow, i);
+			struct short_channel_id_dir scidd;
+
+			get_scidd(rq->gossmap, flow, i, &scidd);
+			struct amount_msat htlc_min =
+			    get_chan_htlc_min(rq, flow->path[i], &scidd);
+			if (amount_msat_less(hop_amt, htlc_min))
+				return false;
+
+			if (!amount_msat_add_fee(&hop_amt, h->base_fee,
+						 h->proportional_fee))
+				abort();
+		}
+	}
+	return true;
+}
+
+static bool check_htlc_max_limits(struct route_query *rq, struct flow **flows)
+{
+
+	for (size_t k = 0; k < tal_count(flows); k++) {
+		struct flow *flow = flows[k];
+		size_t pathlen = tal_count(flow->path);
+		struct amount_msat hop_amt = flow->delivers;
+		for (size_t i = pathlen - 1; i < pathlen; i--) {
+			const struct half_chan *h = flow_edge(flow, i);
+			struct short_channel_id_dir scidd;
+
+			get_scidd(rq->gossmap, flow, i, &scidd);
+			struct amount_msat htlc_max =
+			    get_chan_htlc_max(rq, flow->path[i], &scidd);
+			if (amount_msat_greater(hop_amt, htlc_max))
+				return false;
+
+			if (!amount_msat_add_fee(&hop_amt, h->base_fee,
+						 h->proportional_fee))
+				abort();
+		}
+	}
+	return true;
+}
+
+/* FIXME: add extra constraint maximum route length, use an activation
+ * probability cost for each channel. Recall that every activation cost, eg.
+ * base fee and activation probability can only be properly added modifying the
+ * graph topology by creating an activation node for every half channel. */
+/* FIXME: add extra constraint maximum number of routes, fixes issue 8331. */
+/* FIXME: add a boolean option to make recipient pay for fees, fixes issue 8353.
+ */
+static const char *
+linear_routes(const tal_t *ctx, struct route_query *rq,
+	      struct timemono deadline,
+	      const struct gossmap_node *srcnode,
+	      const struct gossmap_node *dstnode, struct amount_msat amount,
+	      struct amount_msat maxfee, u32 finalcltv, u32 maxdelay,
+	      struct flow ***flows, double *probability,
+	      struct flow **(*solver)(const tal_t *, const struct route_query *,
+				      const struct gossmap_node *,
+				      const struct gossmap_node *,
+				      struct amount_msat, u32, double))
+{
+	const tal_t *working_ctx = tal(ctx, tal_t);
+	const char *error_message;
+	struct amount_msat amount_to_deliver = amount;
+	struct amount_msat feebudget = maxfee;
+
+	/* FIXME: mu is an integer from 0 to MU_MAX that we use to combine fees
+	 * and probability costs, but I think we can make it a real number from
+	 * 0 to 1. */
+	u32 mu = 1;
+	/* we start at 1e-6 and increase it exponentially (x2) up to 10. */
+	double delay_feefactor = 1e-6;
+
+	struct flow **new_flows = NULL;
+	struct amount_msat all_deliver;
+
+	*flows = tal_arr(working_ctx, struct flow *, 0);
+
+	/* Re-use the reservation system to make flows aware of each other. */
+	struct reserve_hop *reservations = new_reservations(working_ctx, rq);
+
+	while (!amount_msat_is_zero(amount_to_deliver)) {
+		if (timemono_after(time_mono(), deadline)) {
+			error_message = rq_log(ctx, rq, LOG_BROKEN,
+					       "%s: timed out after deadline",
+					       __func__);
+			goto fail;
+		}
+
+		new_flows = tal_free(new_flows);
+
+		/* If the amount_to_deliver is very small we better use a single
+		 * path computation because:
+		 * 1. we save cpu cycles
+		 * 2. we have better control over htlc_min violations.
+                 * We need to make the distinction here because after
+                 * refine_with_fees_and_limits we might have a set of flows that
+                 * do not deliver the entire payment amount by just a small
+                 * amount. */
+		if (amount_msat_less_eq(amount_to_deliver,
+					SINGLE_PATH_THRESHOLD)) {
+			new_flows = single_path_flow(working_ctx, rq, srcnode,
+						     dstnode, amount_to_deliver,
+						     mu, delay_feefactor);
+		} else {
+			new_flows =
+			    solver(working_ctx, rq, srcnode, dstnode,
+				   amount_to_deliver, mu, delay_feefactor);
+		}
+
+		if (!new_flows) {
+			error_message = explain_failure(
+			    ctx, rq, srcnode, dstnode, amount_to_deliver);
+			goto fail;
+		}
+
+		error_message =
+			refine_flows(ctx, rq, amount_to_deliver, &new_flows);
+		if (error_message)
+			goto fail;
+
+		/* we finished removing flows and excess */
+		all_deliver = flowset_delivers(rq->plugin, new_flows);
+		if (amount_msat_is_zero(all_deliver)) {
+			/* We removed all flows and we have not modified the
+			 * MCF parameters. We will not have an infinite loop
+			 * here because at least we have disabled some channels.
+			 */
+			continue;
+		}
+
+		/* We might want to overpay sometimes, eg. shadow routing, but
+		 * right now if all_deliver > amount_to_deliver means a bug. */
+		assert(amount_msat_greater_eq(amount_to_deliver, all_deliver));
+
+		/* no flows should send 0 amount */
+		for (size_t i = 0; i < tal_count(new_flows); i++) {
+                        // FIXME: replace all assertions with LOG_BROKEN
+			assert(!amount_msat_is_zero(new_flows[i]->delivers));
+		}
+
+		/* Is this set of flows too expensive?
+		 * We can check if the new flows are within the fee budget,
+		 * however in some cases we have discarded some flows at this
+		 * point and the new flows do not deliver all the value we need
+		 * so that a further solver iteration is needed. Hence we
+		 * check if the fees paid by these new flows are below the
+		 * feebudget proportionally adjusted by the amount this set of
+		 * flows deliver with respect to the total remaining amount,
+		 * ie. we avoid "consuming" all the feebudget if we still need
+		 * to run MCF again for some remaining amount. */
+		struct amount_msat all_fees =
+		    flowset_fee(rq->plugin, new_flows);
+		const double deliver_fraction =
+		    amount_msat_ratio(all_deliver, amount_to_deliver);
+		struct amount_msat partial_feebudget;
+		if (!amount_msat_scale(&partial_feebudget, feebudget,
+				       deliver_fraction)) {
+			error_message =
+			    rq_log(ctx, rq, LOG_BROKEN,
+				   "%s: failed to scale the fee budget (%s) by "
+				   "fraction (%lf)",
+				   __func__, fmt_amount_msat(tmpctx, feebudget),
+				   deliver_fraction);
+			goto fail;
+		}
+		if (amount_msat_greater(all_fees, partial_feebudget)) {
+			if (mu < MU_MAX) {
+				/* all_fees exceed the strong budget limit, try
+				 * to fix it increasing mu. */
+				if (mu == 1)
+					mu = 10;
+				else
+					mu += 10;
+				mu = MIN(mu, MU_MAX);
+				rq_log(
+				    tmpctx, rq, LOG_INFORM,
+				    "The flows had a fee of %s, greater than "
+				    "max of %s, retrying with mu of %u%%...",
+				    fmt_amount_msat(tmpctx, all_fees),
+				    fmt_amount_msat(tmpctx, partial_feebudget),
+				    mu);
+				continue;
+			} else if (amount_msat_greater(all_fees, feebudget)) {
+				/* we cannot increase mu anymore and all_fees
+				 * already exceeds feebudget we fail. */
+				error_message =
+				    rq_log(ctx, rq, LOG_UNUSUAL,
+					   "Could not find route without "
+					   "excessive cost");
+				goto fail;
+			} else {
+				/* mu cannot be increased but at least all_fees
+				 * does not exceed feebudget, we give it a shot.
+				 */
+				rq_log(
+				    tmpctx, rq, LOG_UNUSUAL,
+				    "The flows had a fee of %s, greater than "
+				    "max of %s, but still within the fee "
+				    "budget %s, we accept those flows.",
+				    fmt_amount_msat(tmpctx, all_fees),
+				    fmt_amount_msat(tmpctx, partial_feebudget),
+				    fmt_amount_msat(tmpctx, feebudget));
+			}
+		}
+
+		/* Too much delay? */
+		if (finalcltv + flows_worst_delay(new_flows) > maxdelay) {
+			if (delay_feefactor > 10) {
+				error_message =
+				    rq_log(ctx, rq, LOG_UNUSUAL,
+					   "Could not find route without "
+					   "excessive delays");
+				goto fail;
+			}
+
+			delay_feefactor *= 2;
+			rq_log(tmpctx, rq, LOG_INFORM,
+			       "The worst flow delay is %" PRIu64
+			       " (> %i), retrying with delay_feefactor %f...",
+			       flows_worst_delay(new_flows), maxdelay - finalcltv,
+			       delay_feefactor);
+                        continue;
+		}
+
+		all_fees = AMOUNT_MSAT(0);
+		all_deliver = AMOUNT_MSAT(0);
+		/* add the new flows to the final solution */
+		for (size_t i = 0; i < tal_count(new_flows); i++) {
+			/* last check: every time we add a new reservation to a
+			 * local channel we remove some amount to pay for fees
+			 * on the additional HTLC. */
+			if (create_flow_reservations_verify(rq, &reservations,
+							    new_flows[i])) {
+				tal_arr_expand(flows, new_flows[i]);
+				tal_steal(*flows, new_flows[i]);
+				if (!amount_msat_accumulate(
+					&all_deliver, new_flows[i]->delivers) ||
+				    !amount_msat_accumulate(
+					&all_fees,
+					flow_fee(rq->plugin, new_flows[i])))
+					abort();
+			}
+		}
+
+		if (!amount_msat_deduct(&feebudget, all_fees) ||
+		    !amount_msat_deduct(&amount_to_deliver, all_deliver)) {
+			error_message =
+			    rq_log(ctx, rq, LOG_BROKEN,
+				   "%s: unexpected arithmetic operation "
+				   "failure on amount_msat",
+				   __func__);
+			goto fail;
+		}
+	}
+	/* transfer ownership */
+	*flows = tal_steal(ctx, *flows);
+
+	/* cleanup */
+	tal_free(working_ctx);
+
+	/* all set! Now squash flows that use the same path */
+	squash_flows(ctx, rq, flows);
+
+	/* If we're over the number of parts, try to cram excess into the
+	 * largest-capacity parts */
+	if (tal_count(*flows) > rq->maxparts) {
+		struct amount_msat fee;
+
+		error_message = reduce_num_flows(rq, rq, flows, amount, rq->maxparts);
+		if (error_message) {
+			*flows = tal_free(*flows);
+			return error_message;
+		}
+
+		/* Check fee budget! */
+		fee = flowset_fee(rq->plugin, *flows);
+		if (amount_msat_greater(fee, maxfee)) {
+			error_message = rq_log(rq, rq, LOG_INFORM,
+					       "After reducing the flows to %zu (i.e. maxparts),"
+					       " we had a fee of %s, greater than "
+					       "max of %s.",
+					       tal_count(*flows),
+					       fmt_amount_msat(tmpctx, fee),
+					       fmt_amount_msat(tmpctx, maxfee));
+			if (error_message) {
+				*flows = tal_free(*flows);
+				return error_message;
+			}
+		}
+	}
+
+	/* flows_probability re-does a temporary reservation so we need to call
+	 * it after we have cleaned the reservations we used to build the flows
+	 * hence after we freed working_ctx. */
+	*probability = flows_probability(ctx, rq, flows);
+
+	/* we should have fixed all htlc violations, "don't trust,
+	 * verify" */
+	if (!check_htlc_min_limits(rq, *flows)) {
+		error_message =
+		    rq_log(rq, rq, LOG_BROKEN,
+			   "%s: check_htlc_min_limits failed", __func__);
+		*flows = tal_free(*flows);
+		return error_message;
+	}
+	if (!check_htlc_max_limits(rq, *flows)) {
+		*flows = tal_free(*flows);
+		return rq_log(rq, rq, LOG_BROKEN,
+			      "%s: check_htlc_max_limits failed", __func__);
+	}
+	if (tal_count(*flows) > rq->maxparts) {
+		size_t num_flows = tal_count(*flows);
+		*flows = tal_free(*flows);
+		return rq_log(rq, rq, LOG_BROKEN,
+			      "%s: the number of flows (%zu) exceeds the limit set "
+			      "on payment parts (%" PRIu32
+			      "), please submit a bug report",
+			      __func__, num_flows, rq->maxparts);
+	}
+
+	return NULL;
+fail:
+	/* cleanup */
+	tal_free(working_ctx);
+
+	assert(error_message != NULL);
+	return error_message;
+}
+
+const char *default_routes(const tal_t *ctx, struct route_query *rq,
+			   struct timemono deadline,
+			   const struct gossmap_node *srcnode,
+			   const struct gossmap_node *dstnode,
+			   struct amount_msat amount, struct amount_msat maxfee,
+			   u32 finalcltv, u32 maxdelay, struct flow ***flows,
+			   double *probability)
+{
+	return linear_routes(ctx, rq, deadline, srcnode, dstnode, amount, maxfee,
+			     finalcltv, maxdelay, flows, probability, minflow);
+}
+
+const char *single_path_routes(const tal_t *ctx, struct route_query *rq,
+			       struct timemono deadline,
+			       const struct gossmap_node *srcnode,
+			       const struct gossmap_node *dstnode,
+			       struct amount_msat amount,
+			       struct amount_msat maxfee, u32 finalcltv,
+			       u32 maxdelay, struct flow ***flows,
+			       double *probability)
+{
+	return linear_routes(ctx, rq, deadline, srcnode, dstnode, amount, maxfee,
+			     finalcltv, maxdelay, flows, probability,
+			     single_path_flow);
+}
diff --git a/plugins/askrene/child/mcf.h b/plugins/askrene/child/mcf.h
new file mode 100644
index 00000000..020b72e9
--- /dev/null
+++ b/plugins/askrene/child/mcf.h
@@ -0,0 +1,33 @@
+#ifndef LIGHTNING_PLUGINS_ASKRENE_CHILD_MCF_H
+#define LIGHTNING_PLUGINS_ASKRENE_CHILD_MCF_H
+/* Eduardo Quintela's (lagrang3@protonmail.com) Min Cost Flow implementation
+ * from renepay, as modified to fit askrene */
+#include "config.h"
+#include <ccan/time/time.h>
+#include <common/amount.h>
+#include <common/gossmap.h>
+
+struct route_query;
+
+/* A wrapper to the min. cost flow solver that actually takes into consideration
+ * the extra msats per channel needed to pay for fees. */
+const char *default_routes(const tal_t *ctx, struct route_query *rq,
+			   struct timemono deadline,
+			   const struct gossmap_node *srcnode,
+			   const struct gossmap_node *dstnode,
+			   struct amount_msat amount,
+			   struct amount_msat maxfee, u32 finalcltv,
+			   u32 maxdelay, struct flow ***flows,
+			   double *probability);
+
+/* A wrapper to the single-path constrained solver. */
+const char *single_path_routes(const tal_t *ctx, struct route_query *rq,
+			       struct timemono deadline,
+			       const struct gossmap_node *srcnode,
+			       const struct gossmap_node *dstnode,
+			       struct amount_msat amount,
+			       struct amount_msat maxfee, u32 finalcltv,
+			       u32 maxdelay, struct flow ***flows,
+			       double *probability);
+
+#endif /* LIGHTNING_PLUGINS_ASKRENE_CHILD_MCF_H */
diff --git a/plugins/askrene/child/priorityqueue.c b/plugins/askrene/child/priorityqueue.c
new file mode 100644
index 00000000..3991fb11
--- /dev/null
+++ b/plugins/askrene/child/priorityqueue.c
@@ -0,0 +1,151 @@
+#define NDEBUG 1
+#include "config.h"
+#include <plugins/askrene/child/priorityqueue.h>
+
+/* priorityqueue: a data structure for pairs (key, value) with
+ * 0<=key<max_num_elements, with easy access to elements by key and the pair
+ * with the smallest value. */
+struct priorityqueue {
+	s64 *value;
+	u32 *base;
+	u32 **heapptr;
+	size_t heapsize;
+	struct gheap_ctx gheap_ctx;
+};
+
+static const s64 INFINITE = INT64_MAX;
+
+/* Required a global priorityqueue for gheap. */
+static struct priorityqueue *global_priorityqueue;
+
+/* The heap comparer for priorityqueue search. Since the top element must be the
+ * one with the smallest value, we use the operator >, rather than <. */
+static int priorityqueue_less_comparer(const void *const ctx UNUSED,
+				       const void *const a,
+				       const void *const b) {
+	return global_priorityqueue->value[*(u32 *)a] >
+	       global_priorityqueue->value[*(u32 *)b];
+}
+
+/* The heap move operator for priorityqueue search. */
+static void priorityqueue_item_mover(void *const dst, const void *const src) {
+	u32 src_idx = *(u32 *)src;
+	*(u32 *)dst = src_idx;
+
+	/* we keep track of the pointer position of each element in the heap,
+	 * for easy update. */
+	global_priorityqueue->heapptr[src_idx] = dst;
+}
+
+/* Allocation of resources for the heap. */
+struct priorityqueue *priorityqueue_new(const tal_t *ctx,
+					size_t max_num_nodes) {
+	struct priorityqueue *q = tal(ctx, struct priorityqueue);
+	/* check allocation */
+	if (!q) return NULL;
+
+	q->value = tal_arr(q, s64, max_num_nodes);
+	q->base = tal_arr(q, u32, max_num_nodes);
+	q->heapptr = tal_arrz(q, u32 *, max_num_nodes);
+
+	/* check allocation */
+	if (!q->value || !q->base || !q->heapptr) return tal_free(q);
+
+	q->heapsize = 0;
+	q->gheap_ctx.fanout = 2;
+	q->gheap_ctx.page_chunks = 1024;
+	q->gheap_ctx.item_size = sizeof(q->base[0]);
+	q->gheap_ctx.less_comparer = priorityqueue_less_comparer;
+	q->gheap_ctx.less_comparer_ctx = NULL;
+	q->gheap_ctx.item_mover = priorityqueue_item_mover;
+	return q;
+}
+
+void priorityqueue_init(struct priorityqueue *q) {
+	const size_t max_num_nodes = tal_count(q->value);
+	q->heapsize = 0;
+	for (size_t i = 0; i < max_num_nodes; ++i) {
+		q->value[i] = INFINITE;
+		q->heapptr[i] = NULL;
+	}
+}
+size_t priorityqueue_size(const struct priorityqueue *q) { return q->heapsize; }
+
+size_t priorityqueue_maxsize(const struct priorityqueue *q) {
+	return tal_count(q->value);
+}
+
+static void priorityqueue_append(struct priorityqueue *q, u32 key, s64 value) {
+	assert(priorityqueue_size(q) < priorityqueue_maxsize(q));
+	assert(key < priorityqueue_maxsize(q));
+
+	const size_t pos = q->heapsize;
+
+	q->base[pos] = key;
+	q->value[key] = value;
+	q->heapptr[key] = &(q->base[pos]);
+	q->heapsize++;
+}
+
+void priorityqueue_update(struct priorityqueue *q, u32 key, s64 value) {
+	assert(key < priorityqueue_maxsize(q));
+
+	if (!q->heapptr[key]) {
+		/* not in the heap */
+		priorityqueue_append(q, key, value);
+		global_priorityqueue = q;
+		gheap_restore_heap_after_item_increase(
+		    &q->gheap_ctx, q->base, q->heapsize,
+		    q->heapptr[key] - q->base);
+		global_priorityqueue = NULL;
+		return;
+	}
+
+	if (q->value[key] > value) {
+		/* value decrease */
+		q->value[key] = value;
+
+		global_priorityqueue = q;
+		gheap_restore_heap_after_item_increase(
+		    &q->gheap_ctx, q->base, q->heapsize,
+		    q->heapptr[key] - q->base);
+		global_priorityqueue = NULL;
+	} else {
+		/* value increase */
+		q->value[key] = value;
+
+		global_priorityqueue = q;
+		gheap_restore_heap_after_item_decrease(
+		    &q->gheap_ctx, q->base, q->heapsize,
+		    q->heapptr[key] - q->base);
+		global_priorityqueue = NULL;
+	}
+	/* assert(gheap_is_heap(&q->gheap_ctx,
+	 *                      q->base,
+	 * 		        priorityqueue_size())); */
+}
+
+u32 priorityqueue_top(const struct priorityqueue *q) {
+	assert(!priorityqueue_empty(q));
+	return q->base[0];
+}
+
+bool priorityqueue_empty(const struct priorityqueue *q) {
+	return q->heapsize == 0;
+}
+
+void priorityqueue_pop(struct priorityqueue *q) {
+	if (q->heapsize == 0) return;
+
+	const u32 top = priorityqueue_top(q);
+	assert(q->heapptr[top] == q->base);
+
+	global_priorityqueue = q;
+	gheap_pop_heap(&q->gheap_ctx, q->base, q->heapsize--);
+	global_priorityqueue = NULL;
+	q->heapptr[top] = NULL;
+}
+
+const s64 *priorityqueue_value(const struct priorityqueue *q) {
+	return q->value;
+}
diff --git a/plugins/askrene/child/priorityqueue.h b/plugins/askrene/child/priorityqueue.h
new file mode 100644
index 00000000..93959c86
--- /dev/null
+++ b/plugins/askrene/child/priorityqueue.h
@@ -0,0 +1,35 @@
+#ifndef LIGHTNING_PLUGINS_ASKRENE_CHILD_PRIORITYQUEUE_H
+#define LIGHTNING_PLUGINS_ASKRENE_CHILD_PRIORITYQUEUE_H
+
+/* Defines a priority queue using gheap. */
+
+#include "config.h"
+#include <ccan/short_types/short_types.h>
+#include <ccan/tal/tal.h>
+#include <gheap.h>
+
+/* Allocation of resources for the heap. */
+struct priorityqueue *priorityqueue_new(const tal_t *ctx,
+					size_t max_num_elements);
+
+/* Initialization of the heap for a new priorityqueue search. */
+void priorityqueue_init(struct priorityqueue *priorityqueue);
+
+/* Inserts a new element in the heap. If node_idx was already in the heap then
+ * its value is updated. */
+void priorityqueue_update(struct priorityqueue *priorityqueue, u32 key,
+			  s64 value);
+
+u32 priorityqueue_top(const struct priorityqueue *priorityqueue);
+bool priorityqueue_empty(const struct priorityqueue *priorityqueue);
+void priorityqueue_pop(struct priorityqueue *priorityqueue);
+
+const s64 *priorityqueue_value(const struct priorityqueue *priorityqueue);
+
+/* Number of elements on the heap. */
+size_t priorityqueue_size(const struct priorityqueue *priorityqueue);
+
+/* Maximum number of elements the heap can host */
+size_t priorityqueue_maxsize(const struct priorityqueue *priorityqueue);
+
+#endif /* LIGHTNING_PLUGINS_ASKRENE_CHILD_PRIORITYQUEUE_H */
diff --git a/plugins/askrene/child/refine.c b/plugins/askrene/child/refine.c
new file mode 100644
index 00000000..dd0b0d06
--- /dev/null
+++ b/plugins/askrene/child/refine.c
@@ -0,0 +1,723 @@
+#include "config.h"
+#include <ccan/asort/asort.h>
+#include <ccan/cast/cast.h>
+#include <ccan/tal/str/str.h>
+#include <common/gossmap.h>
+#include <plugins/askrene/askrene.h>
+#include <plugins/askrene/child/flow.h>
+#include <plugins/askrene/child/refine.h>
+#include <plugins/askrene/reserve.h>
+#include <string.h>
+
+/* We (ab)use the reservation system to place temporary reservations
+ * on channels while we are refining each flow.  This has the effect
+ * of making flows aware of each other. */
+
+/* Get the scidd for the i'th hop in flow */
+static void get_scidd(const struct gossmap *gossmap,
+		      const struct flow *flow,
+		      size_t i,
+		      struct short_channel_id_dir *scidd)
+{
+	scidd->scid = gossmap_chan_scid(gossmap, flow->path[i]);
+	scidd->dir = flow->dirs[i];
+}
+
+static void destroy_reservations(struct reserve_hop *rhops, struct askrene *askrene)
+{
+	for (size_t i = 0; i < tal_count(rhops); i++)
+		reserve_remove(askrene->reserved, &rhops[i]);
+}
+
+struct reserve_hop *new_reservations(const tal_t *ctx,
+				     const struct route_query *rq)
+{
+	struct reserve_hop *rhops = tal_arr(ctx, struct reserve_hop, 0);
+
+	/* Unreserve on free */
+	tal_add_destructor2(rhops, destroy_reservations, get_askrene(rq->plugin));
+	return rhops;
+}
+
+static struct reserve_hop *find_reservation(struct reserve_hop *rhops,
+					    const struct short_channel_id_dir *scidd)
+{
+	for (size_t i = 0; i < tal_count(rhops); i++) {
+		if (short_channel_id_dir_eq(scidd, &rhops[i].scidd))
+			return &rhops[i];
+	}
+	return NULL;
+}
+
+/* Add/update reservation: we (ab)use this to temporarily avoid over-usage as
+ * we refine. */
+static void add_reservation(struct reserve_hop **reservations,
+			    const struct route_query *rq,
+			    const struct gossmap_chan *chan,
+			    const struct short_channel_id_dir *scidd,
+			    struct amount_msat amt)
+{
+	struct reserve_hop rhop, *prev;
+	struct askrene *askrene = get_askrene(rq->plugin);
+	size_t idx;
+
+	/* Update in-place if possible */
+	prev = find_reservation(*reservations, scidd);
+	if (prev) {
+		reserve_remove(askrene->reserved, prev);
+		if (!amount_msat_accumulate(&prev->amount, amt))
+			abort();
+		reserve_add(askrene->reserved, prev, rq->cmd->id);
+		return;
+	}
+	rhop.scidd = *scidd;
+	rhop.amount = amt;
+	/* We don't have to restrict it to a layer, since it's transitory:
+	 * nobody else will see this. */
+	rhop.layer = NULL;
+	reserve_add(askrene->reserved, &rhop, rq->cmd->id);
+
+	/* Set capacities entry to 0 so it get_constraints() looks in reserve. */
+	idx = gossmap_chan_idx(rq->gossmap, chan);
+	if (idx < tal_count(rq->capacities))
+		rq->capacities[idx] = 0;
+
+	/* Record so destructor will unreserve */
+	tal_arr_expand(reservations, rhop);
+}
+
+void create_flow_reservations(const struct route_query *rq,
+			      struct reserve_hop **reservations,
+			      const struct flow *flow)
+{
+	struct amount_msat msat;
+
+	msat = flow->delivers;
+	for (int i = tal_count(flow->path) - 1; i >= 0; i--) {
+		const struct half_chan *h = flow_edge(flow, i);
+		struct amount_msat amount_to_reserve;
+		struct short_channel_id_dir scidd;
+
+		get_scidd(rq->gossmap, flow, i, &scidd);
+
+		/* Reserve more for local channels if it reduces capacity */
+		if (!amount_msat_add(&amount_to_reserve, msat,
+				     get_additional_per_htlc_cost(rq, &scidd)))
+			abort();
+
+		add_reservation(reservations, rq, flow->path[i], &scidd,
+				amount_to_reserve);
+		if (!amount_msat_add_fee(&msat,
+					 h->base_fee, h->proportional_fee))
+			plugin_err(rq->plugin, "Adding fee to amount");
+	}
+}
+
+bool create_flow_reservations_verify(const struct route_query *rq,
+				     struct reserve_hop **reservations,
+				     const struct flow *flow)
+{
+	struct amount_msat msat;
+	msat = flow->delivers;
+	for (int i = tal_count(flow->path) - 1; i >= 0; i--) {
+		struct amount_msat known_min, known_max;
+		const struct half_chan *h = flow_edge(flow, i);
+		struct amount_msat amount_to_reserve = msat;
+		struct short_channel_id_dir scidd;
+
+		get_scidd(rq->gossmap, flow, i, &scidd);
+		get_constraints(rq, flow->path[i], flow->dirs[i], &known_min,
+				&known_max);
+		if (amount_msat_greater(amount_to_reserve, known_max))
+			return false;
+
+		if (!amount_msat_add_fee(&msat, h->base_fee,
+					 h->proportional_fee))
+			abort();
+	}
+	create_flow_reservations(rq, reservations, flow);
+	return true;
+}
+
+/* We use an fp16_t approximatin for htlc_max/min: this gets the exact value. */
+static struct amount_msat get_chan_htlc_max(const struct route_query *rq,
+					    const struct gossmap_chan *c,
+					    int dir)
+{
+	struct amount_msat htlc_max;
+
+	gossmap_chan_get_update_details(rq->gossmap,
+					c, dir,
+					NULL, NULL, NULL, NULL, NULL, NULL,
+					NULL, &htlc_max);
+	return htlc_max;
+}
+
+static struct amount_msat get_chan_htlc_min(const struct route_query *rq,
+					    const struct gossmap_chan *c,
+					    int dir)
+{
+	struct amount_msat htlc_min;
+
+	gossmap_chan_get_update_details(rq->gossmap,
+					c, dir,
+					NULL, NULL, NULL, NULL, NULL, NULL,
+					&htlc_min, NULL);
+	return htlc_min;
+}
+
+enum why_capped {
+	CAPPED_HTLC_MAX,
+	CAPPED_CAPACITY,
+};
+
+/* Reverse order: bigger first */
+static int revcmp_flows(struct flow *const *a, struct flow *const *b, void *unused)
+{
+	if (amount_msat_eq((*a)->delivers, (*b)->delivers))
+		return 0;
+	if (amount_msat_greater((*a)->delivers, (*b)->delivers))
+		return -1;
+	return 1;
+}
+
+// TODO: unit test:
+//      -> make a path
+//      -> compute x = flow_max_deliverable
+//      -> check that htlc_max are all satisfied
+//      -> check that (x+1) at least one htlc_max is violated
+/* Given the channel constraints, return the maximum amount that can be
+ * delivered. */
+static struct amount_msat flow_max_deliverable(const struct route_query *rq,
+					       const struct flow *flow)
+{
+	struct amount_msat deliver = AMOUNT_MSAT(-1);
+	for (size_t i = 0; i < tal_count(flow->path); i++) {
+		const struct half_chan *hc = &flow->path[i]->half[flow->dirs[i]];
+		struct amount_msat unused, known_max, htlc_max;
+		deliver = amount_msat_sub_fee(deliver, hc->base_fee,
+					      hc->proportional_fee);
+		htlc_max = get_chan_htlc_max(rq, flow->path[i], flow->dirs[i]);
+		if (amount_msat_greater(deliver, htlc_max))
+			deliver = htlc_max;
+
+		get_constraints(rq, flow->path[i], flow->dirs[i],
+				&unused, &known_max);
+		if (amount_msat_greater(deliver, known_max))
+			deliver = known_max;
+	}
+	return deliver;
+}
+
+// TODO: unit test:
+//      -> make a path
+//      -> compute x = path_min_deliverable
+//      -> check that htlc_min are all satisfied
+//      -> check that (x-1) at least one htlc_min is violated
+/* The least amount that we can deliver at the destination such that when one
+ * computes the hop amounts backwards the htlc_min are always met. */
+static struct amount_msat flow_min_deliverable(const struct route_query *rq,
+					       const struct flow *flow)
+{
+	struct amount_msat least_send = AMOUNT_MSAT(1);
+	const size_t pathlen = tal_count(flow->path);
+
+	for (size_t i = pathlen - 1; i < pathlen; i--) {
+		const struct half_chan *hc = &flow->path[i]->half[flow->dirs[i]];
+		struct amount_msat htlc_min = get_chan_htlc_min(rq, flow->path[i], flow->dirs[i]);
+
+		least_send = amount_msat_max(least_send, htlc_min);
+		if (!amount_msat_add_fee(&least_send, hc->base_fee,
+					 hc->proportional_fee))
+			abort();
+	}
+
+	/* least_send: is the least amount we can send in order to deliver at
+	 * least 1 msat at the destination. */
+	struct amount_msat least_destination = least_send;
+	for (size_t i = 0; i < pathlen; i++) {
+		const struct half_chan *hc = &flow->path[i]->half[flow->dirs[i]];
+		struct amount_msat htlc_min = get_chan_htlc_min(rq, flow->path[i], flow->dirs[i]);
+		struct amount_msat in_value = least_destination;
+		struct amount_msat out_value =
+		    amount_msat_sub_fee(in_value, hc->base_fee,
+					hc->proportional_fee);
+		assert(amount_msat_greater_eq(out_value, htlc_min));
+		struct amount_msat x = out_value;
+		if (!amount_msat_add_fee(&x, hc->base_fee,
+					 hc->proportional_fee))
+			abort();
+		/* if the in_value computed from the out_value is smaller than
+		 * it should, then we add 1msat */
+		if (amount_msat_less(x, in_value) &&
+		    !amount_msat_accumulate(&out_value, AMOUNT_MSAT(1)))
+			abort();
+		/* check conditions */
+		assert(amount_msat_greater_eq(out_value, htlc_min));
+		x = out_value;
+		assert(
+		    amount_msat_add_fee(&x, hc->base_fee,
+					hc->proportional_fee) &&
+		    amount_msat_greater_eq(x, in_value));
+		least_destination = out_value;
+	}
+	return least_destination;
+}
+
+static const char *
+remove_htlc_min_violations(const tal_t *ctx, struct route_query *rq,
+			   const struct flow *flow)
+{
+	const char *error_message = NULL;
+	struct amount_msat msat = flow->delivers;
+	for (size_t i = tal_count(flow->path) - 1; i < tal_count(flow->path);
+	     i--) {
+		struct amount_msat htlc_min = get_chan_htlc_min(rq, flow->path[i], flow->dirs[i]);
+		const struct half_chan *hc = &flow->path[i]->half[flow->dirs[i]];
+		if (amount_msat_less(msat, htlc_min)) {
+			struct short_channel_id_dir scidd;
+			/* FIXME: hoist this! */
+			size_t idx = flow->dirs[i]
+				+ 2 * gossmap_chan_idx(rq->gossmap, flow->path[i]);
+
+			get_scidd(rq->gossmap, flow, i, &scidd);
+			rq_log(
+			    ctx, rq, LOG_INFORM,
+			    "Sending %s across %s would violate htlc_min "
+			    "(~%s), disabling this channel",
+			    fmt_amount_msat(ctx, msat),
+			    fmt_short_channel_id_dir(ctx, &scidd),
+			    fmt_amount_msat(ctx, htlc_min));
+			bitmap_set_bit(rq->disabled_chans, idx);
+			break;
+		}
+		if (!amount_msat_add_fee(
+			&msat, hc->base_fee,
+			hc->proportional_fee)) {
+			error_message =
+			    rq_log(ctx, rq, LOG_BROKEN,
+				   "%s: Adding fee to amount", __func__);
+			break;
+		}
+	}
+	return error_message;
+}
+
+/* Loop over the channels in the path and disable the one with the least
+ * permiting amount based on htlc_max and known max liquidity. */
+static const char *remove_bottleneck(const tal_t *ctx, struct route_query *rq,
+				     const struct flow *flow)
+{
+	const char *error_message = NULL;
+	struct amount_msat min = AMOUNT_MSAT(-1);
+	u32 min_pos = UINT32_MAX;
+	struct amount_msat htlc_max, known_max, unused;
+	struct short_channel_id_dir scidd;
+	size_t idx;
+	for (u32 i = 0; i < tal_count(flow->path); i++) {
+		htlc_max = get_chan_htlc_max(rq, flow->path[i], flow->dirs[i]);
+		get_constraints(rq, flow->path[i], flow->dirs[i], &unused,
+				&known_max);
+		known_max = amount_msat_min(known_max, htlc_max);
+		if (amount_msat_less(known_max, min)) {
+			min = known_max;
+			min_pos = i;
+		}
+	}
+	if (min_pos >= tal_count(flow->path)) {
+		error_message = rq_log(
+		    ctx, rq, LOG_BROKEN,
+		    "%s: failed to find any bottleneck, flow has no hops? %s",
+		    __func__, fmt_flow_full(tmpctx, rq, flow));
+	} else {
+		get_scidd(rq->gossmap, flow, min_pos, &scidd);
+		rq_log(ctx, rq, LOG_INFORM,
+		       "Disabling bottleneck channel %s with "
+		       "htlc_max/known_max at %s",
+		       fmt_short_channel_id_dir(ctx, &scidd),
+		       fmt_amount_msat(ctx, min));
+		idx = flow->dirs[min_pos] +
+		      2 * gossmap_chan_idx(rq->gossmap, flow->path[min_pos]);
+		bitmap_set_bit(rq->disabled_chans, idx);
+	}
+	return error_message;
+}
+
+static struct amount_msat sum_all_deliver(struct flow **flows)
+{
+	struct amount_msat all_deliver = AMOUNT_MSAT(0);
+	for (size_t i = 0; i < tal_count(flows); i++) {
+		if (!amount_msat_accumulate(&all_deliver,
+					    flows[i]->delivers))
+			abort();
+	}
+	return all_deliver;
+}
+
+/* Remove and free the flow */
+static void del_flow_from_arr(struct flow ***flows, size_t i)
+{
+	tal_free((*flows)[i]);
+	tal_arr_remove(flows, i);
+}
+
+/* It reduces the amount of the flows and/or removes some flows in order to
+ * deliver no more than max_deliver. It will leave at least one flow.
+ * Returns the total delivery amount. */
+static struct amount_msat remove_excess(struct flow ***flows,
+					struct amount_msat max_deliver)
+{
+	if (tal_count(*flows) == 0)
+		return AMOUNT_MSAT(0);
+
+	struct amount_msat all_deliver, excess;
+	all_deliver = sum_all_deliver(*flows);
+
+	/* early exit: there is no excess */
+	if (!amount_msat_sub(&excess, all_deliver, max_deliver) ||
+	    amount_msat_is_zero(excess))
+		return all_deliver;
+
+	asort(*flows, tal_count(*flows), revcmp_flows, NULL);
+
+	/* Remove the smaller parts if they deliver less than the
+	 * excess.  */
+	for (int i = tal_count(*flows) - 1; i >= 0; i--) {
+		if (!amount_msat_deduct(&excess,
+				     (*flows)[i]->delivers))
+			break;
+		if (!amount_msat_deduct(&all_deliver,
+				     (*flows)[i]->delivers))
+			abort();
+		del_flow_from_arr(flows, i);
+	}
+
+	/* If we still have some excess, remove it from the
+	 * current flows in the same proportion every flow contributes to the
+	 * total. */
+	struct amount_msat old_excess = excess;
+	struct amount_msat old_deliver = all_deliver;
+	for (size_t i = 0; i < tal_count(*flows); i++) {
+		double fraction = amount_msat_ratio(
+			(*flows)[i]->delivers, old_deliver);
+		struct amount_msat remove;
+
+		if (!amount_msat_scale(&remove, old_excess, fraction))
+			abort();
+
+		/* rounding errors: don't remove more than excess */
+		remove = amount_msat_min(remove, excess);
+
+		if (!amount_msat_deduct(&excess, remove))
+			abort();
+
+		if (!amount_msat_deduct(&all_deliver, remove) ||
+		    !amount_msat_deduct(&(*flows)[i]->delivers, remove))
+			abort();
+	}
+
+	/* any rounding error left, take it from the first */
+	assert(tal_count(*flows) > 0);
+	if (!amount_msat_deduct(&all_deliver, excess) ||
+	    !amount_msat_deduct(&(*flows)[0]->delivers, excess))
+		abort();
+	return all_deliver;
+}
+
+/* Return true (and set shortage) if flow doesn't deliver this much */
+static bool flows_short(struct flow **flows,
+			struct amount_msat deliver,
+			struct amount_msat *shortage)
+{
+	return amount_msat_sub(shortage, deliver, sum_all_deliver(flows))
+		&& !amount_msat_is_zero(*shortage);
+}
+
+/* It increases the flows to meet the deliver target. It does not increase any
+ * flow beyond the tolerance fraction (unless negative).
+ * Returns true if it managed to increase total amount to "deliver". */
+static bool increase_flows(const struct route_query *rq,
+			   struct flow **flows,
+			   struct amount_msat deliver,
+			   double tolerance)
+{
+	const tal_t *working_ctx = tal(NULL, tal_t);
+	struct amount_msat shortage, *ceiling;
+
+	/* Record max we can deliver for each flow, so we don't exceed it */
+	ceiling = tal_arr(working_ctx, struct amount_msat, tal_count(flows));
+	for (size_t i = 0; i < tal_count(flows); i++) {
+		if (tolerance < 0)
+			ceiling[i] = deliver;
+		else if (!amount_msat_scale(&ceiling[i], flows[i]->delivers, 1.0 + tolerance))
+			abort();
+	}
+
+	/* This is naive, but since flows can overlap, increasing one
+	 * can alter the remaining capacity of the others! */
+	while (flows_short(flows, deliver, &shortage)) {
+		size_t best_flownum = 0;
+		struct amount_msat best_remaining = AMOUNT_MSAT(0);
+		struct reserve_hop **reservations;
+		struct amount_msat addition;
+
+		/* Because flows can interact, we reserve them all, removing one at a time. */
+		reservations = tal_arr(NULL, struct reserve_hop *, tal_count(flows));
+		for (size_t i = 0; i < tal_count(flows); i++) {
+			reservations[i] = new_reservations(reservations, rq);
+			create_flow_reservations(rq, &reservations[i], flows[i]);
+		}
+
+		/* Find flow with most excess capacity. */
+		for (size_t i = 0; i < tal_count(flows); i++) {
+			struct amount_msat capacity, remaining;
+
+			/* flow_max_deliverable considers reservations *and*
+			 * htlc_max.  So remove this reservation, to get the
+			 * real maximum for one flow, then replace it. */
+			tal_free(reservations[i]);
+			capacity = flow_max_deliverable(rq, flows[i]);
+			reservations[i] = new_reservations(reservations, rq);
+			create_flow_reservations(rq, &reservations[i], flows[i]);
+
+			/* Don't go above our tolerance */
+			if (amount_msat_greater(capacity, ceiling[i]))
+				capacity = ceiling[i];
+
+			/* We've had a report that this subtract can fail:
+			 * that implies we've pushed a flow past its estimated
+			 * capacity.  That shouldn't happen, but if it does,
+			 * we don't crash */
+			if (!amount_msat_sub(&remaining, capacity, flows[i]->delivers)) {
+			        rq_log(rq, rq, LOG_BROKEN,
+				       "%s: flow %s delivers %s which is more than the path's capacity %s", __func__,
+				       fmt_flow_full(tmpctx, rq, flows[i]),
+				       fmt_amount_msat(tmpctx, flows[i]->delivers),
+				       fmt_amount_msat(tmpctx, capacity));
+				continue;
+			}
+			if (amount_msat_greater(remaining, best_remaining)) {
+				best_flownum = i;
+				best_remaining = remaining;
+			}
+		}
+		tal_free(reservations);
+
+		/* Add 1/n of the remainder, or all we can if that's less than 10 sats. */
+		if (amount_msat_less_sat(shortage, AMOUNT_SAT(10)))
+			addition = shortage;
+		else
+			addition = amount_msat_div_ceil(shortage, tal_count(flows));
+
+		/* Can't add it? */
+		if (amount_msat_less(best_remaining, addition)) {
+			tal_free(working_ctx);
+			return false;
+		}
+
+		if (!amount_msat_accumulate(&flows[best_flownum]->delivers, addition))
+			abort();
+		if (!amount_msat_deduct(&shortage, addition))
+			abort();
+	}
+	tal_free(working_ctx);
+	return true;
+}
+
+const char *refine_flows(const tal_t *ctx, struct route_query *rq,
+			 struct amount_msat deliver, struct flow ***flows)
+{
+	const tal_t *working_ctx = tal(ctx, tal_t);
+	const char *error_message = NULL;
+	struct amount_msat *min_deliverable;
+	size_t *flows_index;
+
+	/* do not deliver more than HTLC_MAX allow us */
+	for (size_t i = 0; i < tal_count(*flows);) {
+		struct amount_msat try_deliver = (*flows)[i]->delivers;
+		struct amount_msat deliverable =
+		    flow_max_deliverable(rq, (*flows)[i]);
+
+		/* We don't expect to have a zero flow amount here. Just report
+		 * it. */
+		if (amount_msat_is_zero(try_deliver)) {
+			rq_log(ctx, rq, LOG_UNUSUAL,
+			       "Tried to refine a flow with zero amount: %s",
+			       fmt_flow_full(tmpctx, rq, (*flows)[i]));
+			del_flow_from_arr(flows, i);
+			continue;
+		}
+
+		/* A path with a very small deliverable amount is not worth the
+		 * effort, and we don't want either a path that for fees and
+		 * HTLC max constraints removes too much from the actual
+		 * delivery amount. In theory the MCF already has partitioned
+		 * the payment in different paths. The refinement step is not
+		 * expected to change the flow by much. */
+		if (amount_msat_less(deliverable, AMOUNT_MSAT(1000)) ||
+		    amount_msat_ratio(deliverable, try_deliver) < 0.2) {
+			error_message = remove_bottleneck(ctx, rq, (*flows)[i]);
+			if (error_message)
+				goto fail;
+			del_flow_from_arr(flows, i);
+			continue;
+		}
+
+		(*flows)[i]->delivers =
+		    amount_msat_min(try_deliver, deliverable);
+		i++;
+	}
+	if (tal_count(*flows) == 0) {
+		/* No flows left to complete the next steps, early exit. */
+		goto fail;
+	}
+
+	/* remove excess from MCF granularity if any */
+	remove_excess(flows, deliver);
+
+	min_deliverable = tal_arrz(working_ctx, struct amount_msat,
+				   tal_count(*flows));
+	flows_index = tal_arrz(working_ctx, size_t, tal_count(*flows));
+	for (size_t i = 0; i < tal_count(*flows); i++) {
+		// FIXME: does flow_max_deliverable work for a single
+		// channel with 0 fees?
+		min_deliverable[i] = flow_min_deliverable(rq, (*flows)[i]);
+		/* We use an array of indexes to keep track of the order
+		 * of the flows. Likewise flows can be removed by simply
+		 * shrinking the flows_index array. */
+		flows_index[i] = i;
+	}
+
+	/* increase flows if necessary to meet the target */
+	increase_flows(rq, *flows, deliver, /* tolerance = */ 0.02);
+
+	/* detect htlc_min violations */
+	for (size_t i = 0; i < tal_count(*flows);) {
+		if (amount_msat_greater_eq((*flows)[i]->delivers,
+					   flow_min_deliverable(rq, (*flows)[i]))) {
+			i++;
+			continue;
+		}
+		error_message = remove_htlc_min_violations(
+		    ctx, rq, (*flows)[i]);
+		if (error_message)
+			goto fail;
+		/* htlc_min is not met for this flow */
+		del_flow_from_arr(flows, i);
+	}
+
+	/* remove 0 amount flows if any */
+	asort(*flows, tal_count(*flows), revcmp_flows, NULL);
+	for (int i = tal_count(*flows) - 1; i >= 0; i--) {
+		if (!amount_msat_is_zero((*flows)[i]->delivers))
+			break;
+		del_flow_from_arr(flows, i);
+	}
+
+	tal_free(working_ctx);
+	return NULL;
+
+fail:
+	tal_free(working_ctx);
+	return error_message;
+}
+
+/* Order of flows in lexicographic order */
+static int cmppath_flows(struct flow *const *a, struct flow *const *b, void *unused)
+{
+	const struct flow *fa = *a, *fb = *b;
+	for (size_t i = 0; i < tal_count(fa->path); i++) {
+		/* Shorter comes first */
+		if (i >= tal_count(fb->path))
+			return 1;
+		if (fa->path[i] < fb->path[i])
+			return -1;
+		if (fa->path[i] > fb->path[i])
+			return 1;
+	}
+	/* fa equal to fb, but is fb longer? */
+	if (tal_count(fb->path) > tal_count(fa->path))
+		return -1;
+	/* equal */
+	return 0;
+}
+
+void squash_flows(const tal_t *ctx, struct route_query *rq,
+		  struct flow ***flows)
+{
+	asort(*flows, tal_count(*flows), cmppath_flows, NULL);
+	for (size_t i = 0; i < tal_count(*flows); i++) {
+		struct flow *flow = (*flows)[i];
+
+		/* same path? We merge */
+		while (i + 1 < tal_count(*flows) &&
+		       cmppath_flows(&flow, &(*flows)[i+1], NULL) == 0) {
+			struct amount_msat combined, max = flow_max_deliverable(rq, flow);
+
+			if (!amount_msat_add(&combined, flow->delivers, (*flows)[i+1]->delivers))
+				abort();
+			/* do we break any HTLC max limits */
+			if (amount_msat_greater(combined, max))
+				break;
+			flow->delivers = combined;
+			del_flow_from_arr(flows, i+1);
+		}
+	}
+}
+
+double flows_probability(const tal_t *ctx, struct route_query *rq,
+			 struct flow ***flows)
+{
+	const tal_t *working_ctx = tal(ctx, tal_t);
+	struct reserve_hop *reservations = new_reservations(working_ctx, rq);
+	double probability = 1.0;
+
+	for (size_t i = 0; i < tal_count(*flows); i++) {
+		probability *= flow_probability((*flows)[i], rq);
+		create_flow_reservations(rq, &reservations, (*flows)[i]);
+	}
+	tal_free(working_ctx);
+	return probability;
+}
+
+const char *reduce_num_flows(const tal_t *ctx,
+			     const struct route_query *rq,
+			     struct flow ***flows,
+			     struct amount_msat deliver,
+			     size_t num_parts)
+{
+	/* Keep the largest flows (not as I originally implemented, the largest
+	 * capacity flows).  Here's Lagrang3's analysis:
+	 *
+	 * I think it is better to keep the largest-deliver flows.  If we only
+	 * go for the highest capacity we may throw away the low cost benefits
+	 * of the MCF.
+
+	 * Hypothetical scenario: MCF finds 3 flows but maxparts=2,
+	 * flow 1: deliver=10, cost=0, capacity=0
+	 * flow 2: deliver=7, cost=1, capacity=5
+	 * flow 3: deliver=1, cost=10, capacity=100
+	 *
+	 * It is better to keep flows 1 and 2 by accomodating 1 more unit of
+	 * flow in flow2 at 1 value expense (per flow), than to keep flows 2 and
+	 * 3 by accomodating 5 more units of flow in flow2 at cost 1 and 5 in
+	 * flow3 at cost 100.
+	 *
+	 * The trade-off is: if we prioritize the delivery value already
+	 * computed by MCF then we find better solutions, but we might fail to
+	 * find feasible solutions sometimes. If we prioritize capacity then we
+	 * generally find bad solutions though we find feasibility more often
+	 * than the alternative.
+	 */
+	size_t orig_num_flows = tal_count(*flows);
+	asort(*flows, orig_num_flows, revcmp_flows, NULL);
+	while (tal_count(*flows) > num_parts)
+		del_flow_from_arr(flows, tal_count(*flows) - 1);
+
+	if (!increase_flows(rq, *flows, deliver, -1.0))
+		return rq_log(ctx, rq, LOG_INFORM,
+			      "Failed to reduce %zu flows down to maxparts (%zu)",
+			      orig_num_flows, num_parts);
+
+	return NULL;
+}
diff --git a/plugins/askrene/child/refine.h b/plugins/askrene/child/refine.h
new file mode 100644
index 00000000..2726e591
--- /dev/null
+++ b/plugins/askrene/child/refine.h
@@ -0,0 +1,44 @@
+#ifndef LIGHTNING_PLUGINS_ASKRENE_CHILD_REFINE_H
+#define LIGHTNING_PLUGINS_ASKRENE_CHILD_REFINE_H
+#include "config.h"
+#include <ccan/tal/tal.h>
+
+struct route_query;
+struct amount_msat;
+struct flow;
+
+struct reserve_hop *new_reservations(const tal_t *ctx,
+				     const struct route_query *rq);
+
+void create_flow_reservations(const struct route_query *rq,
+			      struct reserve_hop **reservations,
+			      const struct flow *flow);
+
+/* create flow reservations, but first verify that the flow indeeds fits in the
+ * liquidity constraints. Takes into account reservations that include per HTLC
+ * extra amounts to pay for onchain fees. */
+bool create_flow_reservations_verify(const struct route_query *rq,
+				     struct reserve_hop **reservations,
+				     const struct flow *flow);
+
+/* Modify flows to meet HTLC min/max requirements.
+ * It takes into account the exact value of the fees expected at each hop.
+ */
+const char *refine_flows(const tal_t *ctx, struct route_query *rq,
+			 struct amount_msat deliver, struct flow ***flows);
+
+/* Duplicated flows are merged into one. This saves in base fee and HTLC fees.
+ */
+void squash_flows(const tal_t *ctx, struct route_query *rq,
+		  struct flow ***flows);
+
+double flows_probability(const tal_t *ctx, struct route_query *rq,
+			 struct flow ***flows);
+
+/* Modify flows so only N remain, if we can.  Returns an error if we cannot. */
+const char *reduce_num_flows(const tal_t *ctx,
+			     const struct route_query *rq,
+			     struct flow ***flows,
+			     struct amount_msat deliver,
+			     size_t num_parts);
+#endif /* LIGHTNING_PLUGINS_ASKRENE_CHILD_REFINE_H */
diff --git a/plugins/askrene/dijkstra.c b/plugins/askrene/dijkstra.c
deleted file mode 100644
index b3f9d39d..00000000
--- a/plugins/askrene/dijkstra.c
+++ /dev/null
@@ -1,186 +0,0 @@
-#define NDEBUG 1
-#include "config.h"
-#include <plugins/askrene/dijkstra.h>
-
-/* In the heap we keep node idx, but in this structure we keep the distance
- * value associated to every node, and their position in the heap as a pointer
- * so that we can update the nodes inside the heap when the distance label is
- * changed.
- *
- * Therefore this is no longer a multipurpose heap, the node_idx must be an
- * index between 0 and less than max_num_nodes. */
-struct dijkstra {
-	//
-	s64 *distance;
-	u32 *base;
-	u32 **heapptr;
-	size_t heapsize;
-	struct gheap_ctx gheap_ctx;
-};
-
-static const s64 INFINITE = INT64_MAX;
-
-/* Required a global dijkstra for gheap. */
-static struct dijkstra *global_dijkstra;
-
-/* The heap comparer for Dijkstra search. Since the top element must be the one
- * with the smallest distance, we use the operator >, rather than <. */
-static int dijkstra_less_comparer(
-	const void *const ctx UNUSED,
-	const void *const a,
-	const void *const b)
-{
-	return global_dijkstra->distance[*(u32*)a]
-		> global_dijkstra->distance[*(u32*)b];
-}
-
-/* The heap move operator for Dijkstra search. */
-static void dijkstra_item_mover(void *const dst, const void *const src)
-{
-	u32 src_idx = *(u32*)src;
-	*(u32*)dst = src_idx;
-
-	// we keep track of the pointer position of each element in the heap,
-	// for easy update.
-	global_dijkstra->heapptr[src_idx] = dst;
-}
-
-/* Allocation of resources for the heap. */
-struct dijkstra *dijkstra_new(const tal_t *ctx, size_t max_num_nodes)
-{
-	struct dijkstra *dijkstra = tal(ctx, struct dijkstra);
-
-	dijkstra->distance = tal_arr(dijkstra,s64,max_num_nodes);
-	dijkstra->base = tal_arr(dijkstra,u32,max_num_nodes);
-	dijkstra->heapptr = tal_arrz(dijkstra,u32*,max_num_nodes);
-
-	dijkstra->heapsize=0;
-
-	dijkstra->gheap_ctx.fanout=2;
-	dijkstra->gheap_ctx.page_chunks=1024;
-	dijkstra->gheap_ctx.item_size=sizeof(dijkstra->base[0]);
-	dijkstra->gheap_ctx.less_comparer=dijkstra_less_comparer;
-	dijkstra->gheap_ctx.less_comparer_ctx=NULL;
-	dijkstra->gheap_ctx.item_mover=dijkstra_item_mover;
-
-	return dijkstra;
-}
-
-
-void dijkstra_init(struct dijkstra *dijkstra)
-{
-	const size_t max_num_nodes = tal_count(dijkstra->distance);
-	dijkstra->heapsize=0;
-	for(size_t i=0;i<max_num_nodes;++i)
-	{
-		dijkstra->distance[i]=INFINITE;
-		dijkstra->heapptr[i] = NULL;
-	}
-}
-size_t dijkstra_size(const struct dijkstra *dijkstra)
-{
-	return dijkstra->heapsize;
-}
-
-size_t dijkstra_maxsize(const struct dijkstra *dijkstra)
-{
-	return tal_count(dijkstra->distance);
-}
-
-static void dijkstra_append(struct dijkstra *dijkstra, u32 node_idx, s64 distance)
-{
-	assert(dijkstra_size(dijkstra) < dijkstra_maxsize(dijkstra));
-	assert(node_idx < dijkstra_maxsize(dijkstra));
-
-	const size_t pos = dijkstra->heapsize;
-
-	dijkstra->base[pos]=node_idx;
-	dijkstra->distance[node_idx]=distance;
-	dijkstra->heapptr[node_idx] = &(dijkstra->base[pos]);
-	dijkstra->heapsize++;
-}
-
-void dijkstra_update(struct dijkstra *dijkstra, u32 node_idx, s64 distance)
-{
-	assert(node_idx < dijkstra_maxsize(dijkstra));
-
-	if(!dijkstra->heapptr[node_idx])
-	{
-		// not in the heap
-		dijkstra_append(dijkstra, node_idx,distance);
-		global_dijkstra = dijkstra;
-		gheap_restore_heap_after_item_increase(
-			&dijkstra->gheap_ctx,
-			dijkstra->base,
-			dijkstra->heapsize,
-			dijkstra->heapptr[node_idx]
-				- dijkstra->base);
-		global_dijkstra = NULL;
-		return;
-	}
-
-	if(dijkstra->distance[node_idx] > distance)
-	{
-		// distance decrease
-		dijkstra->distance[node_idx] = distance;
-
-		global_dijkstra = dijkstra;
-		gheap_restore_heap_after_item_increase(
-			&dijkstra->gheap_ctx,
-			dijkstra->base,
-			dijkstra->heapsize,
-			dijkstra->heapptr[node_idx]
-				- dijkstra->base);
-		global_dijkstra = NULL;
-	}else
-	{
-		// distance increase
-		dijkstra->distance[node_idx] = distance;
-
-		global_dijkstra = dijkstra;
-		gheap_restore_heap_after_item_decrease(
-			&dijkstra->gheap_ctx,
-			dijkstra->base,
-			dijkstra->heapsize,
-			dijkstra->heapptr[node_idx]
-				- dijkstra->base);
-		global_dijkstra = NULL;
-
-	}
-  	// assert(gheap_is_heap(&dijkstra->gheap_ctx,
-	//                      dijkstra->base,
-	// 		     dijkstra_size()));
-}
-
-u32 dijkstra_top(const struct dijkstra *dijkstra)
-{
-	return dijkstra->base[0];
-}
-
-bool dijkstra_empty(const struct dijkstra *dijkstra)
-{
-	return dijkstra->heapsize==0;
-}
-
-void dijkstra_pop(struct dijkstra *dijkstra)
-{
-	if(dijkstra->heapsize==0)
-		return;
-
-	const u32 top = dijkstra_top(dijkstra);
-	assert(dijkstra->heapptr[top]==dijkstra->base);
-
-	global_dijkstra = dijkstra;
-	gheap_pop_heap(
-		&dijkstra->gheap_ctx,
-		dijkstra->base,
-		dijkstra->heapsize--);
-	global_dijkstra = NULL;
-
-	dijkstra->heapptr[top]=NULL;
-}
-
-const s64* dijkstra_distance_data(const struct dijkstra *dijkstra)
-{
-	return dijkstra->distance;
-}
diff --git a/plugins/askrene/dijkstra.h b/plugins/askrene/dijkstra.h
deleted file mode 100644
index f8ff62a8..00000000
--- a/plugins/askrene/dijkstra.h
+++ /dev/null
@@ -1,30 +0,0 @@
-#ifndef LIGHTNING_PLUGINS_ASKRENE_DIJKSTRA_H
-#define LIGHTNING_PLUGINS_ASKRENE_DIJKSTRA_H
-#include "config.h"
-#include <ccan/short_types/short_types.h>
-#include <ccan/tal/tal.h>
-#include <gheap.h>
-
-/* Allocation of resources for the heap. */
-struct dijkstra *dijkstra_new(const tal_t *ctx, size_t max_num_nodes);
-
-/* Initialization of the heap for a new Dijkstra search. */
-void dijkstra_init(struct dijkstra *dijkstra);
-
-/* Inserts a new element in the heap. If node_idx was already in the heap then
- * its distance value is updated. */
-void dijkstra_update(struct dijkstra *dijkstra, u32 node_idx, s64 distance);
-
-u32 dijkstra_top(const struct dijkstra *dijkstra);
-bool dijkstra_empty(const struct dijkstra *dijkstra);
-void dijkstra_pop(struct dijkstra *dijkstra);
-
-const s64* dijkstra_distance_data(const struct dijkstra *dijkstra);
-
-/* Number of elements on the heap. */
-size_t dijkstra_size(const struct dijkstra *dijkstra);
-
-/* Maximum number of elements the heap can host */
-size_t dijkstra_maxsize(const struct dijkstra *dijkstra);
-
-#endif /* LIGHTNING_PLUGINS_ASKRENE_DIJKSTRA_H */
diff --git a/plugins/askrene/explain_failure.c b/plugins/askrene/explain_failure.c
deleted file mode 100644
index be1f137c..00000000
--- a/plugins/askrene/explain_failure.c
+++ /dev/null
@@ -1,324 +0,0 @@
-#include "config.h"
-#include <ccan/tal/str/str.h>
-#include <common/dijkstra.h>
-#include <common/gossmap.h>
-#include <common/route.h>
-#include <plugins/askrene/askrene.h>
-#include <plugins/askrene/explain_failure.h>
-#include <plugins/askrene/layer.h>
-#include <plugins/askrene/reserve.h>
-
-#define NO_USABLE_PATHS_STRING "We could not find a usable set of paths."
-
-/* Dijkstra, reduced to ignore anything but connectivity */
-static bool always_true(const struct gossmap *map,
-			const struct gossmap_chan *c,
-			int dir,
-			struct amount_msat amount,
-			void *unused)
-{
-	return true;
-}
-
-static u64 route_score_one(struct amount_msat fee UNUSED,
-			   struct amount_msat risk UNUSED,
-			   struct amount_msat total UNUSED,
-			   int dir UNUSED,
-			   const struct gossmap_chan *c UNUSED)
-{
-	return 1;
-}
-
-/* This mirrors get_constraints() */
-static const char *why_max_constrained(const tal_t *ctx,
-				       const struct route_query *rq,
-				       struct short_channel_id_dir *scidd,
-				       struct amount_msat amount)
-{
-	char *ret = NULL;
-	const char *reservations;
-	const struct layer *constrains = NULL;
-	struct amount_msat max = amount;
-
-	/* Figure out the layer that constrains us (most) */
-	for (size_t i = 0; i < tal_count(rq->layers); i++) {
-		struct amount_msat min = AMOUNT_MSAT(0), new_max = max;
-
-		layer_apply_constraints(rq->layers[i], scidd, &min, &new_max);
-		if (!amount_msat_eq(new_max, max))
-			constrains = rq->layers[i];
-		max = new_max;
-	}
-
-	if (constrains) {
-		if (!ret)
-			ret = tal_strdup(ctx, "");
-		else
-			tal_append_fmt(&ret, ", ");
-		tal_append_fmt(&ret, "layer %s says max is %s",
-			       layer_name(constrains),
-			       fmt_amount_msat(tmpctx, max));
-	}
-
-	reservations = fmt_reservations(tmpctx, rq->reserved, scidd, rq->layers);
-	if (reservations) {
-		if (!ret)
-			ret = tal_strdup(ctx, "");
-		else
-			tal_append_fmt(&ret, " and ");
-		tal_append_fmt(&ret, "already reserved %s", reservations);
-	}
-
-	/* This seems unlikely, but don't return NULL. */
-	if (!ret)
-		ret = tal_fmt(ctx, "is constrained");
-	return ret;
-}
-
-struct stat {
-	size_t num_channels;
-	struct amount_msat capacity;
-};
-
-struct node_stats {
-	struct stat total, gossip_known, enabled;
-};
-
-enum node_direction {
-	INTO_NODE,
-	OUT_OF_NODE,
-};
-
-static void add_stat(struct stat *stat,
-		     struct amount_msat amount)
-{
-	stat->num_channels++;
-	if (!amount_msat_accumulate(&stat->capacity, amount))
-		abort();
-}
-
-static void node_stats(const struct route_query *rq,
-		       const struct gossmap_node *node,
-		       enum node_direction node_direction,
-		       struct node_stats *stats)
-{
-	memset(stats, 0, sizeof(*stats));
-	for (size_t i = 0; i < node->num_chans; i++) {
-		int dir;
-		struct gossmap_chan *c;
-		struct amount_msat cap_msat;
-
-		c = gossmap_nth_chan(rq->gossmap, node, i, &dir);
-		cap_msat = gossmap_chan_get_capacity(rq->gossmap, c);
-
-		if (node_direction == INTO_NODE)
-			dir = !dir;
-
-		add_stat(&stats->total, cap_msat);
-		if (gossmap_chan_set(c, dir))
-			add_stat(&stats->gossip_known, cap_msat);
-		if (c->half[dir].enabled)
-			add_stat(&stats->enabled, cap_msat);
-	}
-}
-
-static const char *check_capacity(const tal_t *ctx,
-				  const struct route_query *rq,
-				  const struct gossmap_node *node,
-				  enum node_direction node_direction,
-				  struct amount_msat amount,
-				  const char *name)
-{
-	struct node_stats stats;
-
-	node_stats(rq, node, node_direction, &stats);
-	if (amount_msat_greater(amount, stats.total.capacity)) {
-		return rq_log(ctx, rq, LOG_DBG,
-			      NO_USABLE_PATHS_STRING
-			      "  Total %s capacity is only %s"
-			      " (in %zu channels).",
-			      name,
-			      fmt_amount_msat(tmpctx, stats.total.capacity),
-			      stats.total.num_channels);
-	}
-	if (amount_msat_greater(amount, stats.gossip_known.capacity)) {
-		return rq_log(ctx, rq, LOG_DBG,
-			      NO_USABLE_PATHS_STRING
-			      "  Missing gossip for %s: only known %zu/%zu channels, leaving capacity only %s of %s.",
-			      name,
-			      stats.gossip_known.num_channels,
-			      stats.total.num_channels,
-			      fmt_amount_msat(tmpctx, stats.gossip_known.capacity),
-			      fmt_amount_msat(tmpctx, stats.total.capacity));
-	}
-	if (amount_msat_greater(amount, stats.enabled.capacity)) {
-		return rq_log(ctx, rq, LOG_DBG,
-			      NO_USABLE_PATHS_STRING
-			      "  The %s has disabled %zu of %zu channels, leaving capacity only %s of %s.",
-			      name,
-			      stats.total.num_channels - stats.enabled.num_channels,
-			      stats.total.num_channels,
-			      fmt_amount_msat(tmpctx, stats.enabled.capacity),
-			      fmt_amount_msat(tmpctx, stats.total.capacity));
-	}
-	return NULL;
-}
-
-/* Return description of why scidd is disabled scidd */
-static const char *describe_disabled(const tal_t *ctx,
-				     const struct route_query *rq,
-				     const struct gossmap_chan *c,
-				     const struct short_channel_id_dir *scidd)
-{
-	for (int i = tal_count(rq->layers) - 1; i >= 0; i--) {
-		struct gossmap_node *dst = gossmap_nth_node(rq->gossmap, c, !scidd->dir);
-		struct node_id dstid;
-
-		gossmap_node_get_id(rq->gossmap, dst, &dstid);
-		if (layer_disables_node(rq->layers[i], &dstid))
-			return tal_fmt(ctx, "leads to node disabled by layer %s.",
-				       layer_name(rq->layers[i]));
-		else if (layer_disables_chan(rq->layers[i], scidd)) {
-			return tal_fmt(ctx, "marked disabled by layer %s.",
-				       layer_name(rq->layers[i]));
-		}
-	}
-
-	return tal_fmt(ctx, "marked disabled by gossip message.");
-}
-
-static const char *describe_capacity(const tal_t *ctx,
-				     const struct route_query *rq,
-				     const struct short_channel_id_dir *scidd,
-				     struct amount_msat amount)
-{
-	for (int i = tal_count(rq->layers) - 1; i >= 0; i--) {
-		if (layer_created(rq->layers[i], scidd->scid)) {
-			return tal_fmt(ctx, " (created by layer %s) isn't big enough to carry %s.",
-				       layer_name(rq->layers[i]),
-				       fmt_amount_msat(tmpctx, amount));
-		}
-	}
-
-	return tal_fmt(ctx, "isn't big enough to carry %s.",
-		       fmt_amount_msat(tmpctx, amount));
-}
-
-/* We failed to find a flow at all.  Why? */
-const char *explain_failure(const tal_t *ctx,
-			    const struct route_query *rq,
-			    const struct gossmap_node *srcnode,
-			    const struct gossmap_node *dstnode,
-			    struct amount_msat amount)
-{
-	const struct route_hop *hops;
-	const struct dijkstra *dij;
-	char *path;
-	const char *cap_check;
-	const char *explanation;
-	struct short_channel_id_dir scidd;
-	struct gossmap_chan *c;
-	struct amount_msat rolling_amount;
-	struct amount_msat *path_amount;
-
-	/* Do we have enough funds? */
-	cap_check = check_capacity(ctx, rq, srcnode, OUT_OF_NODE,
-				   amount, "source");
-	if (cap_check)
-		return cap_check;
-
-	/* Does destination have enough capacity? */
-	cap_check = check_capacity(ctx, rq, dstnode, INTO_NODE,
-				   amount, "destination");
-	if (cap_check)
-		return cap_check;
-
-	/* OK, fall back to telling them why didn't shortest path
-	 * work.  This covers the "but I have a direct channel!"
-	 * case. */
-	dij = dijkstra(tmpctx, rq->gossmap, dstnode, AMOUNT_MSAT(0), 0,
-		       always_true, route_score_one, NULL);
-	hops = route_from_dijkstra(tmpctx, rq->gossmap, dij, srcnode,
-				   AMOUNT_MSAT(0), 0);
-	if (!hops)
-		return rq_log(ctx, rq, LOG_INFORM,
-			      "There is no connection between source and destination at all");
-
-	/* Description of shortest path */
-	path = tal_strdup(tmpctx, "");
-	for (size_t i = 0; i < tal_count(hops); i++) {
-		tal_append_fmt(&path, "%s%s",
-			       i > 0 ? "->" : "",
-			       fmt_short_channel_id(tmpctx, hops[i].scid));
-	}
-
-	path_amount = tal_arr(tmpctx, struct amount_msat, tal_count(hops));
-	rolling_amount = amount;
-	for (size_t i = tal_count(hops) - 1; i < tal_count(hops); i--) {
-		scidd.scid = hops[i].scid;
-		scidd.dir = hops[i].direction;
-		c = gossmap_find_chan(rq->gossmap, &scidd.scid);
-
-		path_amount[i] = rolling_amount;
-		if (!amount_msat_add_fee(&rolling_amount,
-					 c->half[scidd.dir].base_fee,
-					 c->half[scidd.dir].proportional_fee)) {
-			/* Should not happen, but since the branch exists we use
-			 * it. */
-			explanation = tal_fmt(
-			    tmpctx, "produces a fee overflow for amount %s",
-			    fmt_amount_msat(tmpctx, rolling_amount));
-			return rq_log(ctx, rq, LOG_INFORM,
-				      NO_USABLE_PATHS_STRING
-				      "  The shortest path is %s, but %s %s",
-				      path,
-				      fmt_short_channel_id_dir(tmpctx, &scidd),
-				      explanation);
-		}
-	}
-
-	/* Now walk through this: is it disabled?  Insuff capacity? */
-	for (size_t i = 0; i < tal_count(hops); i++) {
-		struct amount_msat cap_msat, min, max, htlc_max, htlc_min;
-
-		scidd.scid = hops[i].scid;
-		scidd.dir = hops[i].direction;
-		c = gossmap_find_chan(rq->gossmap, &scidd.scid);
-		cap_msat = gossmap_chan_get_capacity(rq->gossmap, c);
-		get_constraints(rq, c, scidd.dir, &min, &max);
-		htlc_max = amount_msat(fp16_to_u64(c->half[scidd.dir].htlc_max));
-		htlc_min = amount_msat(fp16_to_u64(c->half[scidd.dir].htlc_min));
-
-		if (!gossmap_chan_set(c, scidd.dir))
-			explanation = "has no gossip";
-		else if (!c->half[scidd.dir].enabled)
-			explanation = describe_disabled(tmpctx, rq, c, &scidd);
-		else if (amount_msat_greater(path_amount[i], cap_msat))
-			explanation = describe_capacity(tmpctx, rq, &scidd, path_amount[i]);
-		else if (amount_msat_greater(path_amount[i], max))
-			explanation = why_max_constrained(tmpctx, rq,
-							  &scidd, path_amount[i]);
-		else if (amount_msat_greater(path_amount[i], htlc_max))
-			explanation = tal_fmt(tmpctx,
-					      "exceeds htlc_maximum_msat ~%s",
-					      fmt_amount_msat(tmpctx, htlc_max));
-		else if (amount_msat_less(path_amount[i], htlc_min))
-			explanation = tal_fmt(tmpctx,
-					      "below htlc_minumum_msat ~%s",
-					      fmt_amount_msat(tmpctx, htlc_min));
-		else
-			continue;
-
-		return rq_log(ctx, rq, LOG_INFORM,
-			      NO_USABLE_PATHS_STRING
-			      "  The shortest path is %s, but %s %s",
-			      path,
-			      fmt_short_channel_id_dir(tmpctx, &scidd),
-			      explanation);
-	}
-
-	return rq_log(ctx, rq, LOG_BROKEN,
-		       "Actually, I'm not sure why we didn't find the"
-		       " obvious route %s: perhaps this is a bug?",
-		       path);
-}
diff --git a/plugins/askrene/explain_failure.h b/plugins/askrene/explain_failure.h
deleted file mode 100644
index ed470f82..00000000
--- a/plugins/askrene/explain_failure.h
+++ /dev/null
@@ -1,16 +0,0 @@
-#ifndef LIGHTNING_PLUGINS_ASKRENE_EXPLAIN_FAILURE_H
-#define LIGHTNING_PLUGINS_ASKRENE_EXPLAIN_FAILURE_H
-#include "config.h"
-#include <common/amount.h>
-
-struct route_query;
-struct gossmap_node;
-
-/* When MCF returns nothing, try to explain why */
-const char *explain_failure(const tal_t *ctx,
-			    const struct route_query *rq,
-			    const struct gossmap_node *srcnode,
-			    const struct gossmap_node *dstnode,
-			    struct amount_msat amount);
-
-#endif /* LIGHTNING_PLUGINS_ASKRENE_EXPLAIN_FAILURE_H */
diff --git a/plugins/askrene/flow.c b/plugins/askrene/flow.c
deleted file mode 100644
index 9431f7d5..00000000
--- a/plugins/askrene/flow.c
+++ /dev/null
@@ -1,189 +0,0 @@
-#include "config.h"
-#include <assert.h>
-#include <ccan/tal/str/str.h>
-#include <ccan/tal/tal.h>
-#include <common/fp16.h>
-#include <common/overflows.h>
-#include <math.h>
-#include <plugins/askrene/askrene.h>
-#include <plugins/askrene/flow.h>
-#include <plugins/libplugin.h>
-#include <stdio.h>
-
-#ifndef SUPERVERBOSE
-#define SUPERVERBOSE(...)
-#else
-#define SUPERVERBOSE_ENABLED 1
-#endif
-
-/* How much do we deliver to destination using this set of routes */
-struct amount_msat flowset_delivers(struct plugin *plugin,
-				    struct flow **flows)
-{
-	struct amount_msat final = AMOUNT_MSAT(0);
-	for (size_t i = 0; i < tal_count(flows); i++) {
-		if (!amount_msat_accumulate(&final, flows[i]->delivers)) {
-			plugin_err(plugin, "Could not add flowsat %s to %s (%zu/%zu)",
-				   fmt_amount_msat(tmpctx, flows[i]->delivers),
-				   fmt_amount_msat(tmpctx, final),
-				   i, tal_count(flows));
-		}
-	}
-	return final;
-}
-
-/* Stolen whole-cloth from @Lagrang3 in renepay's flow.c.  Wrong
- * because of htlc overhead in reservations! */
-static double edge_probability(const struct route_query *rq,
-			       const struct short_channel_id_dir *scidd,
-			       struct amount_msat sent)
-{
-	struct amount_msat numerator, denominator;
-	struct amount_msat mincap, maxcap, additional;
-	const struct gossmap_chan *c = gossmap_find_chan(rq->gossmap, &scidd->scid);
-
-	get_constraints(rq, c, scidd->dir, &mincap, &maxcap);
-
-	/* We add an extra per-htlc reservation for the *next* HTLC, so we "over-reserve"
-	 * on local channels.  Undo that! */
-	additional = get_additional_per_htlc_cost(rq, scidd);
-	if (!amount_msat_accumulate(&mincap, additional)
-	    || !amount_msat_accumulate(&maxcap, additional))
-		abort();
-
-	if (amount_msat_less_eq(sent, mincap))
-		return 1.0;
-	else if (amount_msat_greater(sent, maxcap))
-		return 0.0;
-
-	/* Linear probability: 1 - (spend - min) / (max - min) */
-
-	/* spend > mincap, from above. */
-	if (!amount_msat_sub(&numerator, sent, mincap))
-		abort();
-	/* This can only fail is maxcap was < mincap,
-	 * so we would be captured above */
-	if (!amount_msat_sub(&denominator, maxcap, mincap))
-		abort();
-	return 1.0 - amount_msat_ratio(numerator, denominator);
-}
-
-struct amount_msat flow_spend(struct plugin *plugin, const struct flow *flow)
-{
-	const size_t pathlen = tal_count(flow->path);
-	struct amount_msat spend = flow->delivers;
-
-	for (int i = (int)pathlen - 1; i >= 0; i--) {
-		const struct half_chan *h = flow_edge(flow, i);
-		if (!amount_msat_add_fee(&spend, h->base_fee,
-					 h->proportional_fee)) {
-			plugin_err(plugin, "Could not add fee %u/%u to amount %s in %i/%zu",
-				   h->base_fee, h->proportional_fee,
-				   fmt_amount_msat(tmpctx, spend),
-				   i, pathlen);
-		}
-	}
-
-	return spend;
-}
-
-struct amount_msat flow_fee(struct plugin *plugin, const struct flow *flow)
-{
-	struct amount_msat spend = flow_spend(plugin, flow);
-	struct amount_msat fee;
-	if (!amount_msat_sub(&fee, spend, flow->delivers)) {
-		plugin_err(plugin, "Could not subtract %s from %s for fee",
-				   fmt_amount_msat(tmpctx, flow->delivers),
-				   fmt_amount_msat(tmpctx, spend));
-	}
-
-	return fee;
-}
-
-struct amount_msat flowset_fee(struct plugin *plugin, struct flow **flows)
-{
-	struct amount_msat fee = AMOUNT_MSAT(0);
-	for (size_t i = 0; i < tal_count(flows); i++) {
-		struct amount_msat this_fee = flow_fee(plugin, flows[i]);
-		if (!amount_msat_accumulate(&fee, this_fee)) {
-			plugin_err(plugin, "Could not add %s to %s for flowset fee",
-				   fmt_amount_msat(tmpctx, this_fee),
-				   fmt_amount_msat(tmpctx, fee));
-		}
-	}
-	return fee;
-}
-
-/* Helper to access the half chan at flow index idx */
-const struct half_chan *flow_edge(const struct flow *flow, size_t idx)
-{
-	assert(flow);
-	assert(idx < tal_count(flow->path));
-	return &flow->path[idx]->half[flow->dirs[idx]];
-}
-
-/* Helper function to find the success_prob for a single flow
- *
- * IMPORTANT: flow->success_prob is misleading, because that's the prob. of
- * success provided that there are no other flows in the current MPP flow set.
- * */
-double flow_probability(const struct flow *flow,
-			const struct route_query *rq)
-{
-	const size_t pathlen = tal_count(flow->path);
-	struct amount_msat spend = flow->delivers;
-	double prob = 1.0;
-
-	for (int i = (int)pathlen - 1; i >= 0; i--) {
-		const struct half_chan *h = flow_edge(flow, i);
-		struct short_channel_id_dir scidd;
-		scidd.scid = gossmap_chan_scid(rq->gossmap, flow->path[i]);
-		scidd.dir = flow->dirs[i];
-
-		prob *= edge_probability(rq, &scidd, spend);
-
-		if (!amount_msat_add_fee(&spend, h->base_fee,
-					 h->proportional_fee)) {
-			plugin_err(rq->plugin, "Could not add fee %u/%u to amount %s in %i/%zu",
-				   h->base_fee, h->proportional_fee,
-				   fmt_amount_msat(tmpctx, spend),
-				   i, pathlen);
-		}
-	}
-
-	return prob;
-}
-
-u64 flow_delay(const struct flow *flow)
-{
-	u64 delay = 0;
-	for (size_t i = 0; i < tal_count(flow->path); i++)
-		delay += flow_edge(flow, i)->delay;
-	return delay;
-}
-
-u64 flows_worst_delay(struct flow **flows)
-{
-	u64 maxdelay = 0;
-	for (size_t i = 0; i < tal_count(flows); i++) {
-		u64 delay = flow_delay(flows[i]);
-		if (delay > maxdelay)
-			maxdelay = delay;
-	}
-	return maxdelay;
-}
-
-const char *fmt_flows_step_scid(const tal_t *ctx,
-				const struct route_query *rq,
-				const struct flow *flow, size_t i)
-{
-	struct short_channel_id_dir scidd;
-
-	scidd.scid = gossmap_chan_scid(rq->gossmap, flow->path[i]);
-	scidd.dir = flow->dirs[i];
-	return fmt_short_channel_id_dir(ctx, &scidd);
-}
-
-#ifndef SUPERVERBOSE_ENABLED
-#undef SUPERVERBOSE
-#endif
diff --git a/plugins/askrene/flow.h b/plugins/askrene/flow.h
deleted file mode 100644
index 28ac7627..00000000
--- a/plugins/askrene/flow.h
+++ /dev/null
@@ -1,69 +0,0 @@
-#ifndef LIGHTNING_PLUGINS_ASKRENE_FLOW_H
-#define LIGHTNING_PLUGINS_ASKRENE_FLOW_H
-#include "config.h"
-#include <bitcoin/short_channel_id.h>
-#include <common/amount.h>
-#include <common/gossmap.h>
-
-struct plugin;
-struct route_query;
-
-/* An actual partial flow. */
-struct flow {
-	const struct gossmap_chan **path;
-	/* The directions to traverse. */
-	int *dirs;
-	/* Amount delivered */
-	struct amount_msat delivers;
-};
-
-/* Helper to access the half chan at flow index idx */
-const struct half_chan *flow_edge(const struct flow *flow, size_t idx);
-
-/* A big number, meaning "don't bother" (not infinite, since you may add) */
-#define FLOW_INF_COST 100000000.0
-
-/* Cost function to send @f msat through @c in direction @dir,
- * given we already have a flow of prev_flow. */
-double flow_edge_cost(const struct gossmap *gossmap,
-		      const struct gossmap_chan *c, int dir,
-		      const struct amount_msat known_min,
-		      const struct amount_msat known_max,
-		      struct amount_msat prev_flow,
-		      struct amount_msat f,
-		      double mu,
-		      double basefee_penalty,
-		      double delay_riskfactor);
-
-/* What's the success probability of this flow in isolation? */
-double flow_probability(const struct flow *flow,
-			const struct route_query *rq);
-
-/* How much do we need to send to make this flow arrive. */
-struct amount_msat flow_spend(struct plugin *plugin, const struct flow *flow);
-
-/* How much do we pay in fees to make this flow arrive. */
-struct amount_msat flow_fee(struct plugin *plugin, const struct flow *flow);
-
-/* What fee to we pay for this entire flow set? */
-struct amount_msat flowset_fee(struct plugin *plugin, struct flow **flows);
-
-/* How much does this entire flowset deliver? */
-struct amount_msat flowset_delivers(struct plugin *plugin,
-				    struct flow **flows);
-
-/* How much CLTV does this flow require? */
-u64 flow_delay(const struct flow *flow);
-
-/* Max CLTV any of these flows requires */
-u64 flows_worst_delay(struct flow **flows);
-
-const char *fmt_flows_step_scid(const tal_t *ctx,
-				const struct route_query *rq,
-				const struct flow *flow, size_t i);
-
-/* When we need to debug */
-const char *fmt_flow_full(const tal_t *ctx,
-			  const struct route_query *rq,
-			  const struct flow *flow);
-#endif /* LIGHTNING_PLUGINS_ASKRENE_FLOW_H */
diff --git a/plugins/askrene/graph.c b/plugins/askrene/graph.c
deleted file mode 100644
index 973e6adc..00000000
--- a/plugins/askrene/graph.c
+++ /dev/null
@@ -1,65 +0,0 @@
-#include "config.h"
-#include <plugins/askrene/graph.h>
-
-/* in the background add the actual arc or dual arc */
-static void graph_push_outbound_arc(struct graph *graph, const struct arc arc,
-				    const struct node node)
-{
-	assert(arc.idx < graph_max_num_arcs(graph));
-	assert(node.idx < graph_max_num_nodes(graph));
-
-	/* arc is already added, skip */
-	if (graph->arc_tail[arc.idx].idx != INVALID_INDEX)
-		return;
-
-	graph->arc_tail[arc.idx] = node;
-
-	const struct arc first_arc = graph->node_adjacency_first[node.idx];
-	graph->node_adjacency_next[arc.idx] = first_arc;
-	graph->node_adjacency_first[node.idx] = arc;
-}
-
-bool graph_add_arc(struct graph *graph, const struct arc arc,
-		   const struct node from, const struct node to)
-{
-	assert(from.idx < graph->max_num_nodes);
-	assert(to.idx < graph->max_num_nodes);
-
-	const struct arc dual = arc_dual(graph, arc);
-
-	if (arc.idx >= graph->max_num_arcs || dual.idx >= graph->max_num_arcs)
-		return false;
-
-	graph_push_outbound_arc(graph, arc, from);
-	graph_push_outbound_arc(graph, dual, to);
-
-	return true;
-}
-
-struct graph *graph_new(const tal_t *ctx, const size_t max_num_nodes,
-			const size_t max_num_arcs, const size_t arc_dual_bit)
-{
-	struct graph *graph;
-	graph = tal(ctx, struct graph);
-
-	graph->max_num_arcs = max_num_arcs;
-	graph->max_num_nodes = max_num_nodes;
-	graph->arc_dual_bit = arc_dual_bit;
-
-	graph->arc_tail = tal_arr(graph, struct node, graph->max_num_arcs);
-	graph->node_adjacency_first =
-	    tal_arr(graph, struct arc, graph->max_num_nodes);
-	graph->node_adjacency_next =
-	    tal_arr(graph, struct arc, graph->max_num_arcs);
-
-	/* initialize with invalid indexes so that we know these slots have
-	 * never been used, eg. arc/node is newly created */
-	for (size_t i = 0; i < graph->max_num_arcs; i++)
-		graph->arc_tail[i] = node_obj(INVALID_INDEX);
-	for (size_t i = 0; i < graph->max_num_nodes; i++)
-		graph->node_adjacency_first[i] = arc_obj(INVALID_INDEX);
-	for (size_t i = 0; i < graph->max_num_nodes; i++)
-		graph->node_adjacency_next[i] = arc_obj(INVALID_INDEX);
-
-	return graph;
-}
diff --git a/plugins/askrene/graph.h b/plugins/askrene/graph.h
deleted file mode 100644
index e84ed131..00000000
--- a/plugins/askrene/graph.h
+++ /dev/null
@@ -1,171 +0,0 @@
-#ifndef LIGHTNING_PLUGINS_ASKRENE_GRAPH_H
-#define LIGHTNING_PLUGINS_ASKRENE_GRAPH_H
-
-/* Defines a graph data structure. */
-
-#include "config.h"
-#include <assert.h>
-#include <ccan/short_types/short_types.h>
-#include <ccan/tal/tal.h>
-
-#define INVALID_INDEX 0xffffffff
-
-/* A directed arc in a graph.
- * It is a simple data object for typesafey. */
-struct arc {
-	/* arc's index */
-	u32 idx;
-};
-
-/* A node in a graph.
- * It is a simple data object for typesafety. */
-struct node {
-	/* node's index */
-	u32 idx;
-};
-
-static inline struct arc arc_obj(u32 index)
-{
-	struct arc arc = {.idx = index};
-	return arc;
-}
-static inline struct node node_obj(u32 index)
-{
-	struct node node = {.idx = index};
-	return node;
-}
-
-/* A graph's topology. */
-struct graph {
-	/* Every arc emanates from a node, the tail.
-	 * The head of the arc is the tail of the dual. */
-	struct node *arc_tail;
-
-	/* Adjacency data for nodes. Used to move in a graph in the direction of
-	 * the arcs by looping over all arcs that exit a node.
-	 *
-	 * For every directed arc there is a dual in the opposite direction,
-	 * therefore we can use the same adjacency information to traverse in
-	 * the head to tails direction as well. */
-	struct arc *node_adjacency_next;
-	struct arc *node_adjacency_first;
-
-	size_t max_num_arcs, max_num_nodes;
-
-	/* Bit that must be flipped to obtain the dual of an arc. */
-	size_t arc_dual_bit;
-};
-
-//////////////////////////////////////////////////////////////////////////////
-
-static inline size_t graph_max_num_arcs(const struct graph *graph)
-{
-	return graph->max_num_arcs;
-}
-static inline size_t graph_max_num_nodes(const struct graph *graph)
-{
-	return graph->max_num_nodes;
-}
-
-/* Give me the dual of an arc. */
-static inline struct arc arc_dual(const struct graph *graph, struct arc arc)
-{
-	arc.idx ^= (1U << graph->arc_dual_bit);
-	return arc;
-}
-
-/* Is this arc a dual? */
-static inline bool arc_is_dual(const struct graph *graph, struct arc arc)
-{
-	return (arc.idx & (1U << graph->arc_dual_bit)) != 0;
-}
-
-/* Give me the node at the tail of an arc. */
-static inline struct node arc_tail(const struct graph *graph,
-				   const struct arc arc)
-{
-	assert(arc.idx < graph_max_num_arcs(graph));
-	return graph->arc_tail[arc.idx];
-}
-
-/* Give me the node at the head of an arc. */
-static inline struct node arc_head(const struct graph *graph,
-				   const struct arc arc)
-{
-	const struct arc dual = arc_dual(graph, arc);
-	assert(dual.idx < graph_max_num_arcs(graph));
-	return graph->arc_tail[dual.idx];
-}
-
-/* We use an arc array but not all arcs in that array do exist in the graph. */
-static inline bool arc_enabled(const struct graph *graph, const struct arc arc)
-{
-	return graph->arc_tail[arc.idx].idx < graph->max_num_nodes;
-}
-
-/* Used to loop over the arcs that exit a node.
- *
- * for example:
- *
- * void show(struct graph *graph, struct node node) {
- * 	printf("Showing node %" PRIu32 "\n", node.idx);
- * 	for (struct arc arc = node_adjacency_begin(graph, node);
- * 	     !node_adjacency_end(arc);
- * 	     arc = node_adjacency_next(graph, arc)) {
- * 		printf("arc id: %" PRIu32 ", (%" PRIu32 " -> %" PRIu32 ")\n",
- * 		       arc.idx,
- * 		       arc_tail(graph, arc).idx,
- * 		       arc_head(graph, arc).idx);
- * 	}
- * }
- * */
-static inline struct arc node_adjacency_begin(const struct graph *graph,
-					      const struct node node)
-{
-	assert(node.idx < graph_max_num_nodes(graph));
-	return graph->node_adjacency_first[node.idx];
-}
-static inline bool node_adjacency_end(const struct arc arc)
-{
-	return arc.idx == INVALID_INDEX;
-}
-static inline struct arc node_adjacency_next(const struct graph *graph,
-					     const struct arc arc)
-{
-	assert(arc.idx < graph_max_num_arcs(graph));
-	return graph->node_adjacency_next[arc.idx];
-}
-
-/* Used to loop over the arcs that enter a node. */
-static inline struct arc node_rev_adjacency_begin(const struct graph *graph,
-						  const struct node node)
-{
-	return arc_dual(graph, node_adjacency_begin(graph, node));
-}
-static inline bool node_rev_adjacency_end(const struct arc arc)
-{
-	return arc.idx == INVALID_INDEX;
-}
-static inline struct arc node_rev_adjacency_next(const struct graph *graph,
-						 const struct arc arc)
-{
-	return arc_dual(graph,
-			node_adjacency_next(graph, arc_dual(graph, arc)));
-}
-
-/* This call adds an arc to the graph, it adds also the dual automatically.
- * An arc cannot be added twice, if the caller tries to do add the same arc
- * twice the second call is ignored.
- * The call fails if the arc or its dual do not fit into max_num_arcs. */
-bool graph_add_arc(struct graph *graph, const struct arc arc,
-		   const struct node from, const struct node to);
-
-/* Creates a graph object. Nodes and arcs are indexed from 0 to max_num_nodes-1
- * and max_num_arcs-1 respectively. The max_num_arcs should be big enough to
- * accomodate also the dual arcs, ie. if the maximum index for a problem arc is
- * I then Idual = I^(1<<arc_dual_bit) must be a valid arc index
- * Idual<max_num_arcs. */
-struct graph *graph_new(const tal_t *ctx, const size_t max_num_nodes,
-			const size_t max_num_arcs, const size_t arc_dual_bit);
-
-#endif /* LIGHTNING_PLUGINS_ASKRENE_GRAPH_H */
diff --git a/plugins/askrene/layer.h b/plugins/askrene/layer.h
index 55931cb9..ab9ab71e 100644
--- a/plugins/askrene/layer.h
+++ b/plugins/askrene/layer.h
@@ -14,6 +14,7 @@
 #include <common/node_id.h>
 
 struct askrene;
+struct command;
 struct layer;
 struct json_stream;
 
diff --git a/plugins/askrene/mcf.c b/plugins/askrene/mcf.c
deleted file mode 100644
index c3bf28fc..00000000
--- a/plugins/askrene/mcf.c
+++ /dev/null
@@ -1,1644 +0,0 @@
-#include "config.h"
-#include <assert.h>
-#include <ccan/asort/asort.h>
-#include <ccan/bitmap/bitmap.h>
-#include <ccan/list/list.h>
-#include <ccan/tal/str/str.h>
-#include <ccan/tal/tal.h>
-#include <common/utils.h>
-#include <float.h>
-#include <inttypes.h>
-#include <math.h>
-#include <plugins/askrene/algorithm.h>
-#include <plugins/askrene/askrene.h>
-#include <plugins/askrene/dijkstra.h>
-#include <plugins/askrene/explain_failure.h>
-#include <plugins/askrene/flow.h>
-#include <plugins/askrene/graph.h>
-#include <plugins/askrene/mcf.h>
-#include <plugins/askrene/refine.h>
-#include <plugins/libplugin.h>
-#include <stdint.h>
-
-/* # Optimal payments
- *
- * In this module we reduce the routing optimization problem to a linear
- * cost optimization problem and find a solution using MCF algorithms.
- * The optimization of the routing itself doesn't need a precise numerical
- * solution, since we can be happy near optimal results; e.g. paying 100 msat or
- * 101 msat for fees doesn't make any difference if we wish to deliver 1M sats.
- * On the other hand, we are now also considering Pickhard's
- * [1] model to improve payment reliability,
- * hence our optimization moves to a 2D space: either we like to maximize the
- * probability of success of a payment or minimize the routing fees, or
- * alternatively we construct a function of the two that gives a good compromise.
- *
- * Therefore from now own, the definition of optimal is a matter of choice.
- * To simplify the API of this module, we think the best way to state the
- * problem is:
- *
- * 	Find a routing solution that pays the least of fees while keeping
- * 	the probability of success above a certain value `min_probability`.
- *
- *
- * # Fee Cost
- *
- * Routing fees is non-linear function of the payment flow x, that's true even
- * without the base fee:
- *
- * 	fee_msat = base_msat + floor(millionths*x_msat / 10^6)
- *
- * We approximate this fee into a linear function by computing a slope `c_fee` such
- * that:
- *
- * 	fee_microsat = c_fee * x_sat
- *
- * Function `linear_fee_cost` computes `c_fee` based on the base and
- * proportional fees of a channel.
- * The final product if microsat because if only
- * the proportional fee was considered we can have c_fee = millionths.
- * Moving to costs based in msats means we have to either truncate payments
- * below 1ksats or estimate as 0 cost for channels with less than 1000ppm.
- *
- * TODO(eduardo): shall we build a linear cost function in msats?
- *
- * # Probability cost
- *
- * The probability of success P of the payment is the product of the prob. of
- * success of forwarding parts of the payment over all routing channels. This
- * problem is separable if we log it, and since we would like to increase P,
- * then we can seek to minimize -log(P), and that's our prob. cost function [1].
- *
- * 	- log P = sum_{i} - log P_i
- *
- * The probability of success `P_i` of sending some flow `x` on a channel with
- * liquidity l in the range a<=l<b is
- *
- * 	P_{a,b}(x) = (b-x)/(b-a); for x > a
- * 		   = 1.         ; for x <= a
- *
- * Notice that unlike the similar formula in [1], the one we propose does not
- * contain the quantization shot noise for counting states. The formula remains
- * valid independently of the liquidity units (sats or msats).
- *
- * The cost associated to probability P is then -k log P, where k is some
- * constant. For k=1 we get the following table:
- *
- * 	prob | cost
- * 	-----------
- * 	0.01 | 4.6
- * 	0.02 | 3.9
- * 	0.05 | 3.0
- * 	0.10 | 2.3
- * 	0.20 | 1.6
- * 	0.50 | 0.69
- * 	0.80 | 0.22
- * 	0.90 | 0.10
- * 	0.95 | 0.05
- * 	0.98 | 0.02
- * 	0.99 | 0.01
- *
- * Clearly -log P(x) is non-linear; we try to linearize it piecewise:
- * split the channel into 4 arcs representing 4 liquidity regions:
- *
- * 	arc_0 -> [0, a)
- * 	arc_1 -> [a, a+(b-a)*f1)
- * 	arc_2 -> [a+(b-a)*f1, a+(b-a)*f2)
- * 	arc_3 -> [a+(b-a)*f2, a+(b-a)*f3)
- *
- * where f1 = 0.5, f2 = 0.8, f3 = 0.95;
- * We fill arc_0's capacity with complete certainty P=1, then if more flow is
- * needed we start filling the capacity in arc_1 until the total probability
- * of success reaches P=0.5, then arc_2 until P=1-0.8=0.2, and finally arc_3 until
- * P=1-0.95=0.05. We don't go further than 5% prob. of success per channel.
-
- * TODO(eduardo): this channel linearization is hard coded into
- * `CHANNEL_PIVOTS`, maybe we can parametrize this to take values from the config file.
- *
- * With this choice, the slope of the linear cost function becomes:
- *
- * 	m_0 = 0
- * 	m_1 = 1.38 k /(b-a)
- * 	m_2 = 3.05 k /(b-a)
- * 	m_3 = 9.24 k /(b-a)
- *
- * Notice that one of the assumptions in [2] for the MCF problem is that flows
- * and the slope of the costs functions are integer numbers. The only way we
- * have at hand to make it so, is to choose a universal value of `k` that scales
- * up the slopes so that floor(m_i) is not zero for every arc.
- *
- * # Combine fee and prob. costs
- *
- * We attempt to solve the original problem of finding the solution that
- * pays the least fees while keeping the prob. of success above a certain value,
- * by constructing a cost function which is a linear combination of fee and
- * prob. costs.
- * TODO(eduardo): investigate how this procedure is justified,
- * possibly with the use of Lagrange optimization theory.
- *
- * At first, prob. and fee costs live in different dimensions, they cannot be
- * summed, it's like comparing apples and oranges.
- * However we propose to scale the prob. cost by a global factor k that
- * translates into the monetization of prob. cost.
- *
- * This was chosen empirically from examination of typical network values.
- *
- * # References
- *
- * [1] Pickhardt and Richter, https://arxiv.org/abs/2107.05322
- * [2] R.K. Ahuja, T.L. Magnanti, and J.B. Orlin. Network Flows:
- * Theory, Algorithms, and Applications. Prentice Hall, 1993.
- *
- *
- * TODO(eduardo) it would be interesting to see:
- * how much do we pay for reliability?
- * Cost_fee(most reliable solution) - Cost_fee(cheapest solution)
- *
- * TODO(eduardo): it would be interesting to see:
- * how likely is the most reliable path with respect to the cheapest?
- * Prob(reliable)/Prob(cheapest) = Exp(Cost_prob(cheapest)-Cost_prob(reliable))
- *
- * */
-
-#define PARTS_BITS 2
-#define CHANNEL_PARTS (1 << PARTS_BITS)
-
-// These are the probability intervals we use to decompose a channel into linear
-// cost function arcs.
-static const double CHANNEL_PIVOTS[]={0,0.5,0.8,0.95};
-
-static const s64 INFINITE = INT64_MAX;
-static const s64 MU_MAX = 100;
-
-/* every payment under 1000sat will be routed through a single path */
-static const struct amount_msat SINGLE_PATH_THRESHOLD = AMOUNT_MSAT(1000000);
-
-/* Let's try this encoding of arcs:
- * Each channel `c` has two possible directions identified by a bit
- * `half` or `!half`, and each one of them has to be
- * decomposed into 4 liquidity parts in order to
- * linearize the cost function, but also to solve MCF
- * problem we need to keep track of flows in the
- * residual network hence we need for each directed arc
- * in the network there must be another arc in the
- * opposite direction refered to as it's dual. In total
- * 1+2+1 additional bits of information:
- *
- * 	(chan_idx)(half)(part)(dual)
- *
- * That means, for each channel we need to store the
- * information of 16 arcs. If we implement a convex-cost
- * solver then we can reduce that number to size(half)size(dual)=4.
- *
- * In the adjacency of a `node` we are going to store
- * the outgoing arcs. If we ever need to loop over the
- * incoming arcs then we will define a reverse adjacency
- * API.
- * Then for each outgoing channel `(c,half)` there will
- * be 4 parts for the actual residual capacity, hence
- * with the dual bit set to 0:
- *
- * 	(c,half,0,0)
- * 	(c,half,1,0)
- * 	(c,half,2,0)
- * 	(c,half,3,0)
- *
- * and also we need to consider the dual arcs
- * corresponding to the channel direction `(c,!half)`
- * (the dual has reverse direction):
- *
- * 	(c,!half,0,1)
- * 	(c,!half,1,1)
- * 	(c,!half,2,1)
- * 	(c,!half,3,1)
- *
- * These are the 8 outgoing arcs relative to `node` and
- * associated with channel `c`. The incoming arcs will
- * be:
- *
- * 	(c,!half,0,0)
- * 	(c,!half,1,0)
- * 	(c,!half,2,0)
- * 	(c,!half,3,0)
- *
- * 	(c,half,0,1)
- * 	(c,half,1,1)
- * 	(c,half,2,1)
- * 	(c,half,3,1)
- *
- * but they will be stored as outgoing arcs on the peer
- * node `next`.
- *
- * I hope this will clarify my future self when I forget.
- *
- * */
-
-/*
- * We want to use the whole number here for convenience, but
- * we can't us a union, since bit order is implementation-defined and
- * we want chanidx on the highest bits:
- *
- * [ 0       1 2     3            4 5 6 ... 31 ]
- *   dual    part    chandir      chanidx
- */
-#define ARC_DUAL_BITOFF (0)
-#define ARC_PART_BITOFF (1)
-#define ARC_CHANDIR_BITOFF (1 + PARTS_BITS)
-#define ARC_CHANIDX_BITOFF (1 + PARTS_BITS + 1)
-#define ARC_CHANIDX_BITS  (32 - ARC_CHANIDX_BITOFF)
-
-/* How many arcs can we have for a single channel?
- * linearization parts, both directions, and dual */
-#define ARCS_PER_CHANNEL ((size_t)1 << (PARTS_BITS + 1 + 1))
-
-static inline void arc_to_parts(struct arc arc,
-				u32 *chanidx,
-				int *chandir,
-				u32 *part,
-				bool *dual)
-{
-	if (chanidx)
-		*chanidx = (arc.idx >> ARC_CHANIDX_BITOFF);
-	if (chandir)
-		*chandir = (arc.idx >> ARC_CHANDIR_BITOFF) & 1;
-	if (part)
-		*part = (arc.idx >> ARC_PART_BITOFF) & ((1 << PARTS_BITS)-1);
-	if (dual)
-		*dual = (arc.idx >> ARC_DUAL_BITOFF) & 1;
-}
-
-static inline struct arc arc_from_parts(u32 chanidx, int chandir, u32 part, bool dual)
-{
-	struct arc arc;
-
-	assert(part < CHANNEL_PARTS);
-	assert(chandir == 0 || chandir == 1);
-	assert(chanidx < (1U << ARC_CHANIDX_BITS));
-	arc.idx = ((u32)dual << ARC_DUAL_BITOFF)
-		| (part << ARC_PART_BITOFF)
-		| ((u32)chandir << ARC_CHANDIR_BITOFF)
-		| (chanidx << ARC_CHANIDX_BITOFF);
-	return arc;
-}
-
-#define MAX(x, y) (((x) > (y)) ? (x) : (y))
-#define MIN(x, y) (((x) < (y)) ? (x) : (y))
-
-struct pay_parameters {
-	const struct route_query *rq;
-	const struct gossmap_node *source;
-	const struct gossmap_node *target;
-
-	// how much we pay
-	struct amount_msat amount;
-
-	/* base unit for computation, ie. accuracy */
-	struct amount_msat accuracy;
-
-	// channel linearization parameters
-	double cap_fraction[CHANNEL_PARTS],
-	       cost_fraction[CHANNEL_PARTS];
-
-	double delay_feefactor;
-	double base_fee_penalty;
-};
-
-/* Helper function.
- * Given an arc of the network (not residual) give me the flow. */
-static s64 get_arc_flow(
-		const s64 *arc_residual_capacity,
-		const struct graph *graph,
-		const struct arc arc)
-{
-	assert(!arc_is_dual(graph, arc));
-	struct arc dual = arc_dual(graph, arc);
-	assert(dual.idx < tal_count(arc_residual_capacity));
-	return arc_residual_capacity[dual.idx];
-}
-
-/* Set *capacity to value, up to *cap_on_capacity.  Reduce cap_on_capacity */
-static void set_capacity(s64 *capacity, u64 value, u64 *cap_on_capacity)
-{
-	*capacity = MIN(value, *cap_on_capacity);
-	*cap_on_capacity -= *capacity;
-}
-
-/* Helper to check whether a channel is available */
-static bool channel_is_available(const struct route_query *rq,
-				 const struct gossmap_chan *chan, const int dir)
-{
-	const u32 c_idx = gossmap_chan_idx(rq->gossmap, chan);
-	return gossmap_chan_set(chan, dir) && chan->half[dir].enabled &&
-	       !bitmap_test_bit(rq->disabled_chans, c_idx * 2 + dir);
-}
-
-/* FIXME: unit test this */
-/* The probability of forwarding a payment amount given a high and low liquidity
- * bounds.
- * @low: the liquidity is known to be greater or equal than "low"
- * @high: the liquidity is known to be less than "high"
- * @amount: how much is required to forward */
-static double pickhardt_richter_probability(struct amount_msat low,
-					    struct amount_msat high,
-					    struct amount_msat amount)
-{
-	struct amount_msat all_states, good_states;
-	if (amount_msat_greater_eq(amount, high))
-		return 0.0;
-	if (!amount_msat_deduct(&amount, low))
-		return 1.0;
-	if (!amount_msat_sub(&all_states, high, low))
-		abort(); // we expect high > low
-	if (!amount_msat_sub(&good_states, all_states, amount))
-		abort(); // we expect high > amount
-	return amount_msat_ratio(good_states, all_states);
-}
-
-// TODO(eduardo): unit test this
-/* Split a directed channel into parts with linear cost function. */
-static void linearize_channel(const struct pay_parameters *params,
-			      const struct gossmap_chan *c, const int dir,
-			      s64 *capacity, double *cost)
-{
-	struct amount_msat mincap, maxcap;
-
-	/* This takes into account any payments in progress. */
-	get_constraints(params->rq, c, dir, &mincap, &maxcap);
-
-	/* Assume if min > max, min is wrong */
-	if (amount_msat_greater(mincap, maxcap))
-		mincap = maxcap;
-
-	u64 a = amount_msat_ratio_floor(mincap, params->accuracy),
-	    b = 1 + amount_msat_ratio_floor(maxcap, params->accuracy);
-
-	/* An extra bound on capacity, here we use it to reduce the flow such
-	 * that it does not exceed htlcmax.
-	 * The cap con capacity is not greater than the amount of payment units
-	 * (msat/accuracy). The way a channel is decomposed into linear cost
-	 * arcs (code below) in ascending cost order ensures that the only the
-	 * necessary capacity to forward the payment is allocated in the lower
-	 * cost arcs. This may lead to some arcs in the decomposition (at the
-	 * high cost end) to have a capacity of 0, and we can prune them while
-	 * keeping the solution optimal. */
-	u64 cap_on_capacity =
-	    MIN(amount_msat_ratio_floor(gossmap_chan_htlc_max(c, dir),
-					params->accuracy),
-		amount_msat_ratio_ceil(params->amount, params->accuracy));
-
-	set_capacity(&capacity[0], a, &cap_on_capacity);
-	cost[0]=0;
-	for(size_t i=1;i<CHANNEL_PARTS;++i)
-	{
-		set_capacity(&capacity[i], params->cap_fraction[i]*(b-a), &cap_on_capacity);
-
-		cost[i] = params->cost_fraction[i] * 1000
-			  * amount_msat_ratio(params->amount, params->accuracy)
-			  / (b - a);
-	}
-}
-
-static int cmp_u64(const u64 *a, const u64 *b, void *unused)
-{
-	if (*a < *b)
-		return -1;
-	if (*a > *b)
-		return 1;
-	return 0;
-}
-
-static int cmp_double(const double *a, const double *b, void *unused)
-{
-	if (*a < *b)
-		return -1;
-	if (*a > *b)
-		return 1;
-	return 0;
-}
-
-static double get_median_ratio(const tal_t *working_ctx,
-			       const struct graph *graph,
-			       const double *arc_prob_cost,
-			       const s64 *arc_fee_cost)
-{
-	const size_t max_num_arcs = graph_max_num_arcs(graph);
-	u64 *u64_arr = tal_arr(working_ctx, u64, max_num_arcs);
-	double *double_arr = tal_arr(working_ctx, double, max_num_arcs);
-	size_t n = 0;
-
-	for (struct arc arc = {.idx=0};arc.idx < max_num_arcs; ++arc.idx) {
-		/* scan real arcs, not unused id slots or dual arcs */
-		if (arc_is_dual(graph, arc) || !arc_enabled(graph, arc))
-			continue;
-		assert(n < max_num_arcs/2);
-		u64_arr[n] = arc_fee_cost[arc.idx];
-		double_arr[n] = arc_prob_cost[arc.idx];
-		n++;
-	}
-	asort(u64_arr, n, cmp_u64, NULL);
-	asort(double_arr, n, cmp_double, NULL);
-
-	/* Empty network, or tiny probability, nobody cares */
-	if (n == 0 || double_arr[n/2] < 0.001)
-		return 1;
-
-	/* You need to scale arc_prob_cost by this to match arc_fee_cost */
-	return u64_arr[n/2] / double_arr[n/2];
-}
-
-static void combine_cost_function(const tal_t *working_ctx,
-				  const struct graph *graph,
-				  const double *arc_prob_cost,
-				  const s64 *arc_fee_cost, const s8 *biases,
-				  s64 mu, s64 *arc_cost)
-{
-	/* probabilty and fee costs are not directly comparable!
-	 * Scale by ratio of (positive) medians. */
-	const double k =
-	    get_median_ratio(working_ctx, graph, arc_prob_cost, arc_fee_cost);
-	const double ln_30 = log(30);
-	const size_t max_num_arcs = graph_max_num_arcs(graph);
-
-	for(struct arc arc = {.idx=0};arc.idx < max_num_arcs; ++arc.idx)
-	{
-		if (arc_is_dual(graph, arc) || !arc_enabled(graph, arc))
-			continue;
-
-		const double pcost = arc_prob_cost[arc.idx];
-		const s64 fcost = arc_fee_cost[arc.idx];
-		double combined;
-		u32 chanidx;
-		int chandir;
-		s32 bias;
-
-		assert(fcost != INFINITE);
-		assert(pcost != DBL_MAX);
-		combined = fcost*mu + (MU_MAX-mu)*pcost*k;
-
-		/* Bias is in human scale, where "bigger is better" */
-		arc_to_parts(arc, &chanidx, &chandir, NULL, NULL);
-		bias = biases[(chanidx << 1) | chandir];
-		if (bias != 0) {
-			/* After some trial and error, this gives a nice
-			 * dynamic range (25 seems to be "infinite" in
-			 * practice):
-			 *    e^(-bias / (100/ln(30)))
-			 */
-			double bias_factor = exp(-bias / (100 / ln_30));
-			arc_cost[arc.idx] = combined * bias_factor;
-		} else {
-			arc_cost[arc.idx] = combined;
-		}
-		/* and the respective dual */
-		struct arc dual = arc_dual(graph, arc);
-		arc_cost[dual.idx] = -combined;
-	}
-}
-
-/* Get the fee cost associated to this directed channel.
- * Cost is expressed as PPM of the payment.
- *
- * Choose and integer `c_fee` to linearize the following fee function
- *
- *  	fee_msat = base_msat + floor(millionths*x_msat / 10^6)
- *
- * into
- *
- *  	fee = c_fee/10^6 * x
- *
- *  use `base_fee_penalty` to weight the base fee and `delay_feefactor` to
- *  weight the CLTV delay.
- *  */
-static s64 linear_fee_cost(u32 base_fee, u32 proportional_fee, u16 cltv_delta,
-			   double base_fee_penalty,
-			   double delay_feefactor)
-{
-	s64 pfee = proportional_fee,
-	    bfee = base_fee,
-	    delay = cltv_delta;
-
-	return pfee + bfee* base_fee_penalty+ delay*delay_feefactor;
-}
-
-/* This is inversely proportional to the amount we expect to send.  Let's
- * assume we will send ~10th of the total amount per path.  But note
- * that it converts to parts per million! */
-static double base_fee_penalty_estimate(struct amount_msat amount)
-{
-	return amount_msat_ratio(AMOUNT_MSAT(10000000), amount);
-}
-
-static void init_linear_network(const tal_t *ctx,
-				const struct pay_parameters *params,
-				struct graph **graph, double **arc_prob_cost,
-				s64 **arc_fee_cost, s64 **arc_capacity)
-{
-	const struct gossmap *gossmap = params->rq->gossmap;
-	const size_t max_num_chans = gossmap_max_chan_idx(gossmap);
-	const size_t max_num_arcs = max_num_chans * ARCS_PER_CHANNEL;
-	const size_t max_num_nodes = gossmap_max_node_idx(gossmap);
-
-	*graph = graph_new(ctx, max_num_nodes, max_num_arcs, ARC_DUAL_BITOFF);
-	*arc_prob_cost = tal_arr(ctx, double, max_num_arcs);
-	for (size_t i = 0; i < max_num_arcs; ++i)
-		(*arc_prob_cost)[i] = DBL_MAX;
-
-	*arc_fee_cost = tal_arr(ctx, s64, max_num_arcs);
-	for (size_t i = 0; i < max_num_arcs; ++i)
-		(*arc_fee_cost)[i] = INT64_MAX;
-
-	*arc_capacity = tal_arrz(ctx, s64, max_num_arcs);
-
-	for(struct gossmap_node *node = gossmap_first_node(gossmap);
-	    node;
-	    node=gossmap_next_node(gossmap,node))
-	{
-		const u32 node_id = gossmap_node_idx(gossmap,node);
-
-		for(size_t j=0;j<node->num_chans;++j)
-		{
-			int half;
-			const struct gossmap_chan *c = gossmap_nth_chan(gossmap,
-			                                                node, j, &half);
-
-			if (!channel_is_available(params->rq, c, half))
-				continue;
-
-			/* If a channel insists on more than our total, remove it */
-			if (amount_msat_less(params->amount, gossmap_chan_htlc_min(c, half)))
-				continue;
-
-			const u32 chan_id = gossmap_chan_idx(gossmap, c);
-
-			const struct gossmap_node *next = gossmap_nth_node(gossmap,
-									   c,!half);
-
-			const u32 next_id = gossmap_node_idx(gossmap,next);
-
-			if(node_id==next_id)
-				continue;
-
-			// `cost` is the word normally used to denote cost per
-			// unit of flow in the context of MCF.
-			double prob_cost[CHANNEL_PARTS];
-			s64 capacity[CHANNEL_PARTS];
-
-			// split this channel direction to obtain the arcs
-			// that are outgoing to `node`
-			linearize_channel(params, c, half, capacity, prob_cost);
-
-			/* linear fee_cost per unit of flow */
-			const s64 fee_cost = linear_fee_cost(
-				c->half[half].base_fee,
-				c->half[half].proportional_fee,
-				c->half[half].delay,
-				params->base_fee_penalty,
-				params->delay_feefactor);
-
-			// let's subscribe the 4 parts of the channel direction
-			// (c,half), the dual of these guys will be subscribed
-			// when the `i` hits the `next` node.
-			for(size_t k=0;k<CHANNEL_PARTS;++k)
-			{
-				/* prune arcs with 0 capacity */
-				if (capacity[k] == 0)
-					continue;
-
-				struct arc arc = arc_from_parts(chan_id, half, k, false);
-
-				graph_add_arc(*graph, arc,
-					      node_obj(node_id),
-					      node_obj(next_id));
-
-				(*arc_capacity)[arc.idx] = capacity[k];
-				(*arc_prob_cost)[arc.idx] = prob_cost[k];
-				(*arc_fee_cost)[arc.idx] = fee_cost;
-
-				// + the respective dual
-				struct arc dual = arc_dual(*graph, arc);
-
-				(*arc_capacity)[dual.idx] = 0;
-				(*arc_prob_cost)[dual.idx] = -prob_cost[k];
-				(*arc_fee_cost)[dual.idx] = -fee_cost;
-			}
-		}
-	}
-}
-
-// flow on directed channels
-struct chan_flow
-{
-	s64 half[2];
-};
-
-/* Search in the network a path of positive flow until we reach a node with
- * positive balance (returns a node idx with positive balance)
- * or we discover a cycle (returns a node idx with 0 balance).
- * */
-static struct node find_path_or_cycle(
-		const tal_t *working_ctx,
-		const struct route_query *rq,
-		const struct chan_flow *chan_flow,
-		const struct node source,
-		const s64 *balance,
-
-		const struct gossmap_chan **prev_chan,
-		int *prev_dir,
-		u32 *prev_idx)
-{
-	const struct gossmap *gossmap = rq->gossmap;
-	const size_t max_num_nodes = gossmap_max_node_idx(gossmap);
-	bitmap *visited =
-	    tal_arrz(working_ctx, bitmap, BITMAP_NWORDS(max_num_nodes));
-	u32 final_idx = source.idx;
-	bitmap_set_bit(visited, final_idx);
-
-	/* It is guaranteed to halt, because we either find a node with
-	 * balance[]>0 or we hit a node twice and we stop. */
-	while (balance[final_idx] <= 0) {
-		u32 updated_idx = INVALID_INDEX;
-		struct gossmap_node *cur =
-		    gossmap_node_byidx(gossmap, final_idx);
-
-		for (size_t i = 0; i < cur->num_chans; ++i) {
-			int dir;
-			const struct gossmap_chan *c =
-			    gossmap_nth_chan(gossmap, cur, i, &dir);
-
-			if (!channel_is_available(rq, c, dir))
-				continue;
-
-			const u32 c_idx = gossmap_chan_idx(gossmap, c);
-
-			/* follow the flow */
-			if (chan_flow[c_idx].half[dir] > 0) {
-				const struct gossmap_node *n =
-				    gossmap_nth_node(gossmap, c, !dir);
-				u32 next_idx = gossmap_node_idx(gossmap, n);
-
-				prev_dir[next_idx] = dir;
-				prev_chan[next_idx] = c;
-				prev_idx[next_idx] = final_idx;
-
-				updated_idx = next_idx;
-				break;
-			}
-		}
-
-		assert(updated_idx != INVALID_INDEX);
-		assert(updated_idx != final_idx);
-		final_idx = updated_idx;
-
-		if (bitmap_test_bit(visited, updated_idx)) {
-			/* We have seen this node before, we've found a cycle.
-			 */
-			assert(balance[updated_idx] <= 0);
-			break;
-		}
-		bitmap_set_bit(visited, updated_idx);
-	}
-	return node_obj(final_idx);
-}
-
-struct list_data
-{
-	struct list_node list;
-	struct flow *flow_path;
-};
-
-/* Given a path from a node with negative balance to a node with positive
- * balance, compute the bigest flow and substract it from the nodes balance and
- * the channels allocation. */
-static struct flow *substract_flow(const tal_t *ctx,
-				   const struct pay_parameters *params,
-				   const struct node source,
-				   const struct node sink,
-				   s64 *balance, struct chan_flow *chan_flow,
-				   const u32 *prev_idx, const int *prev_dir,
-				   const struct gossmap_chan *const *prev_chan)
-{
-	const struct gossmap *gossmap = params->rq->gossmap;
-	assert(balance[source.idx] < 0);
-	assert(balance[sink.idx] > 0);
-	s64 delta = -balance[source.idx];
-	size_t length = 0;
-	delta = MIN(delta, balance[sink.idx]);
-
-	/* We can only walk backwards, now get me the legth of the path and the
-	 * max flow we can send through this route. */
-	for (u32 cur_idx = sink.idx; cur_idx != source.idx;
-	     cur_idx = prev_idx[cur_idx]) {
-		assert(cur_idx != INVALID_INDEX);
-		const int dir = prev_dir[cur_idx];
-		const struct gossmap_chan *const chan = prev_chan[cur_idx];
-
-		/* we could optimize here by caching the idx of the channels in
-		 * the path, but the bottleneck of the algorithm is the MCF
-		 * computation not here. */
-		const u32 chan_idx = gossmap_chan_idx(gossmap, chan);
-
-		delta = MIN(delta, chan_flow[chan_idx].half[dir]);
-		length++;
-	}
-
-	struct flow *f = tal(ctx, struct flow);
-	f->path = tal_arr(f, const struct gossmap_chan *, length);
-	f->dirs = tal_arr(f, int, length);
-
-	/* Walk again and substract the flow value (delta). */
-	assert(delta > 0);
-	balance[source.idx] += delta;
-	balance[sink.idx] -= delta;
-	for (u32 cur_idx = sink.idx; cur_idx != source.idx;
-	     cur_idx = prev_idx[cur_idx]) {
-		const int dir = prev_dir[cur_idx];
-		const struct gossmap_chan *const chan = prev_chan[cur_idx];
-		const u32 chan_idx = gossmap_chan_idx(gossmap, chan);
-
-		length--;
-		/* f->path and f->dirs contain the channels in the path in the
-		 * correct order. */
-		f->path[length] = chan;
-		f->dirs[length] = dir;
-
-		chan_flow[chan_idx].half[dir] -= delta;
-	}
-	if (!amount_msat_mul(&f->delivers, params->accuracy, delta))
-		abort();
-	return f;
-}
-
-/* Substract a flow cycle from the channel allocation. */
-static void substract_cycle(const struct gossmap *gossmap,
-			    const struct node sink,
-			    struct chan_flow *chan_flow, const u32 *prev_idx,
-			    const int *prev_dir,
-			    const struct gossmap_chan *const *prev_chan)
-{
-	s64 delta = INFINITE;
-	u32 cur_idx;
-
-	/* Compute greatest flow in this cycle. */
-	for (cur_idx = sink.idx; cur_idx!=INVALID_INDEX;) {
-		const int dir = prev_dir[cur_idx];
-		const struct gossmap_chan *const chan = prev_chan[cur_idx];
-		const u32 chan_idx = gossmap_chan_idx(gossmap, chan);
-
-		delta = MIN(delta, chan_flow[chan_idx].half[dir]);
-
-		cur_idx = prev_idx[cur_idx];
-		if (cur_idx == sink.idx)
-			/* we have come back full circle */
-			break;
-	}
-	assert(cur_idx==sink.idx);
-
-	/* Walk again and substract the flow value (delta). */
-	assert(delta < INFINITE);
-	assert(delta > 0);
-
-	for (cur_idx = sink.idx;cur_idx!=INVALID_INDEX;) {
-		const int dir = prev_dir[cur_idx];
-		const struct gossmap_chan *const chan = prev_chan[cur_idx];
-		const u32 chan_idx = gossmap_chan_idx(gossmap, chan);
-
-		chan_flow[chan_idx].half[dir] -= delta;
-
-		cur_idx = prev_idx[cur_idx];
-		if (cur_idx == sink.idx)
-			/* we have come back full circle */
-			break;
-	}
-	assert(cur_idx==sink.idx);
-}
-
-/* Given a flow in the residual network, build a set of payment flows in the
- * gossmap that corresponds to this flow. */
-static struct flow **
-get_flow_paths(const tal_t *ctx,
-	       const tal_t *working_ctx,
-	       const struct pay_parameters *params,
-	       const struct graph *graph,
-	       const s64 *arc_residual_capacity)
-{
-	struct flow **flows = tal_arr(ctx,struct flow*,0);
-
-	const size_t max_num_chans = gossmap_max_chan_idx(params->rq->gossmap);
-	struct chan_flow *chan_flow = tal_arrz(working_ctx,struct chan_flow,max_num_chans);
-
-	const size_t max_num_nodes = gossmap_max_node_idx(params->rq->gossmap);
-	s64 *balance = tal_arrz(working_ctx,s64,max_num_nodes);
-
-	const struct gossmap_chan **prev_chan
-		= tal_arr(working_ctx,const struct gossmap_chan *,max_num_nodes);
-
-
-	int *prev_dir = tal_arr(working_ctx,int,max_num_nodes);
-	u32 *prev_idx = tal_arr(working_ctx, u32, max_num_nodes);
-
-	for (u32 node_idx = 0; node_idx < max_num_nodes; node_idx++)
-		prev_idx[node_idx] = INVALID_INDEX;
-
-	// Convert the arc based residual network flow into a flow in the
-	// directed channel network.
-	// Compute balance on the nodes.
-	for (struct node n = {.idx = 0}; n.idx < max_num_nodes; n.idx++) {
-		for(struct arc arc = node_adjacency_begin(graph,n);
-		        !node_adjacency_end(arc);
-			arc = node_adjacency_next(graph,arc))
-		{
-			if(arc_is_dual(graph, arc))
-				continue;
-			struct node m = arc_head(graph,arc);
-			s64 flow = get_arc_flow(arc_residual_capacity,
-						graph, arc);
-			u32 chanidx;
-			int chandir;
-
-			balance[n.idx] -= flow;
-			balance[m.idx] += flow;
-
-			arc_to_parts(arc, &chanidx, &chandir, NULL, NULL);
-			chan_flow[chanidx].half[chandir] +=flow;
-		}
-	}
-
-	// Select all nodes with negative balance and find a flow that reaches a
-	// positive balance node.
-	for (struct node source = {.idx = 0}; source.idx < max_num_nodes;
-	     source.idx++) {
-		// this node has negative balance, flows leaves from here
-		while (balance[source.idx] < 0) {
-			prev_chan[source.idx] = NULL;
-			struct node sink = find_path_or_cycle(
-			    working_ctx, params->rq, chan_flow, source,
-			    balance, prev_chan, prev_dir, prev_idx);
-
-			if (balance[sink.idx] > 0)
-			/* case 1. found a path */
-			{
-				struct flow *fp = substract_flow(
-				    flows, params, source, sink, balance,
-				    chan_flow, prev_idx, prev_dir, prev_chan);
-
-				tal_arr_expand(&flows, fp);
-			} else
-			/* case 2. found a cycle */
-			{
-				substract_cycle(params->rq->gossmap, sink, chan_flow,
-						prev_idx, prev_dir, prev_chan);
-			}
-		}
-	}
-	return flows;
-}
-
-/* Given a single path build a flow set. */
-static struct flow **
-get_flow_singlepath(const tal_t *ctx, const struct pay_parameters *params,
-		    const struct graph *graph, const struct gossmap *gossmap,
-		    const struct node source, const struct node destination,
-		    const u64 pay_amount, const struct arc *prev)
-{
-	struct flow **flows, *f;
-	flows = tal_arr(ctx, struct flow *, 1);
-	f = flows[0] = tal(flows, struct flow);
-
-	size_t length = 0;
-
-	for (u32 cur_idx = destination.idx; cur_idx != source.idx;) {
-		assert(cur_idx != INVALID_INDEX);
-		length++;
-		struct arc arc = prev[cur_idx];
-		struct node next = arc_tail(graph, arc);
-		cur_idx = next.idx;
-	}
-	f->path = tal_arr(f, const struct gossmap_chan *, length);
-	f->dirs = tal_arr(f, int, length);
-
-	for (u32 cur_idx = destination.idx; cur_idx != source.idx;) {
-		int chandir;
-		u32 chanidx;
-		struct arc arc = prev[cur_idx];
-		arc_to_parts(arc, &chanidx, &chandir, NULL, NULL);
-
-		length--;
-		f->path[length] = gossmap_chan_byidx(gossmap, chanidx);
-		f->dirs[length] = chandir;
-
-		struct node next = arc_tail(graph, arc);
-		cur_idx = next.idx;
-	}
-	f->delivers = params->amount;
-	return flows;
-}
-
-// TODO(eduardo): choose some default values for the minflow parameters
-/* eduardo: I think it should be clear that this module deals with linear
- * flows, ie. base fees are not considered. Hence a flow along a path is
- * described with a sequence of directed channels and one amount.
- * In the `pay_flow` module there are dedicated routes to compute the actual
- * amount to be forward on each hop.
- *
- * TODO(eduardo): notice that we don't pay fees to forward payments with local
- * channels and we can tell with absolute certainty the liquidity on them.
- * Check that local channels have fee costs = 0 and bounds with certainty (min=max). */
-// TODO(eduardo): we should LOG_DBG the process of finding the MCF while
-// adjusting the frugality factor.
-static struct flow **minflow(const tal_t *ctx,
-			     const struct route_query *rq,
-			     const struct gossmap_node *source,
-			     const struct gossmap_node *target,
-			     struct amount_msat amount,
-			     u32 mu,
-			     double delay_feefactor)
-{
-	struct flow **flow_paths;
-	/* We allocate everything off this, and free it at the end,
-	 * as we can be called multiple times without cleaning tmpctx! */
-	tal_t *working_ctx = tal(NULL, char);
-	struct pay_parameters *params = tal(working_ctx, struct pay_parameters);
-
-	params->rq = rq;
-	params->source = source;
-	params->target = target;
-	params->amount = amount;
-	/* -> We reduce the granularity of the flow by limiting the subdivision
-	 * of the payment amount into 1000 units of flow. That reduces the
-	 * computational burden for algorithms that depend on it, eg. "capacity
-	 * scaling" and "successive shortest path".
-	 * -> Using Ceil operation instead of Floor so that
-	 *      accuracy x 1000 >= amount
-	 * */
-	params->accuracy =
-	    amount_msat_max(AMOUNT_MSAT(1), amount_msat_div_ceil(amount, 1000));
-
-	// template the channel partition into linear arcs
-	params->cap_fraction[0]=0;
-	params->cost_fraction[0]=0;
-	for(size_t i =1;i<CHANNEL_PARTS;++i)
-	{
-		params->cap_fraction[i]=CHANNEL_PIVOTS[i]-CHANNEL_PIVOTS[i-1];
-		params->cost_fraction[i]=
-			log((1-CHANNEL_PIVOTS[i-1])/(1-CHANNEL_PIVOTS[i]))
-			/params->cap_fraction[i];
-	}
-
-	params->delay_feefactor = delay_feefactor;
-	params->base_fee_penalty = base_fee_penalty_estimate(amount);
-
-	// build the uncertainty network with linearization and residual arcs
-	struct graph *graph;
-	double *arc_prob_cost;
-	s64 *arc_fee_cost;
-	s64 *arc_capacity;
-	init_linear_network(working_ctx, params, &graph, &arc_prob_cost,
-			    &arc_fee_cost, &arc_capacity);
-
-	const size_t max_num_arcs = graph_max_num_arcs(graph);
-	const size_t max_num_nodes = graph_max_num_nodes(graph);
-	s64 *arc_cost;
-	s64 *node_potential;
-	s64 *node_excess;
-	arc_cost = tal_arrz(working_ctx, s64, max_num_arcs);
-	node_potential = tal_arrz(working_ctx, s64, max_num_nodes);
-	node_excess = tal_arrz(working_ctx, s64, max_num_nodes);
-
-	const struct node dst = {.idx = gossmap_node_idx(rq->gossmap, target)};
-	const struct node src = {.idx = gossmap_node_idx(rq->gossmap, source)};
-
-
-	/* Since we have constraint accuracy, ask to find a payment solution
-	 * that can pay a bit more than the actual value rathen than undershoot it.
-	 * That's why we use the ceil function here. */
-	const u64 pay_amount =
-	    amount_msat_ratio_ceil(params->amount, params->accuracy);
-
-	if (!simple_feasibleflow(working_ctx, graph, src, dst,
-				 arc_capacity, pay_amount)) {
-		rq_log(tmpctx, rq, LOG_INFORM,
-		       "%s failed: unable to find a feasible flow.", __func__);
-		goto fail;
-	}
-	combine_cost_function(working_ctx, graph, arc_prob_cost, arc_fee_cost,
-			      rq->biases, mu, arc_cost);
-
-	/* We solve a linear MCF problem. */
-	if (!mcf_refinement(working_ctx,
-			    graph,
-			    node_excess,
-			    arc_capacity,
-			    arc_cost,
-			    node_potential)) {
-		rq_log(tmpctx, rq, LOG_BROKEN,
-		       "%s: MCF optimization step failed", __func__);
-		goto fail;
-	}
-
-	/* We dissect the solution of the MCF into payment routes.
-	 * Actual amounts considering fees are computed for every
-	 * channel in the routes. */
-	flow_paths = get_flow_paths(ctx, working_ctx, params,
-				    graph, arc_capacity);
-	if(!flow_paths){
-		rq_log(tmpctx, rq, LOG_BROKEN,
-		       "%s: failed to extract flow paths from the MCF solution",
-		       __func__);
-		goto fail;
-	}
-	tal_free(working_ctx);
-	return flow_paths;
-
-fail:
-	tal_free(working_ctx);
-	return NULL;
-}
-
-/* Initialize the data vectors for the single-path solver. */
-static void init_linear_network_single_path(
-    const tal_t *ctx, const struct pay_parameters *params, struct graph **graph,
-    double **arc_prob_cost, s64 **arc_fee_cost, s64 **arc_capacity)
-{
-	const size_t max_num_chans = gossmap_max_chan_idx(params->rq->gossmap);
-	const size_t max_num_arcs = max_num_chans * ARCS_PER_CHANNEL;
-	const size_t max_num_nodes = gossmap_max_node_idx(params->rq->gossmap);
-
-	*graph = graph_new(ctx, max_num_nodes, max_num_arcs, ARC_DUAL_BITOFF);
-	*arc_prob_cost = tal_arr(ctx, double, max_num_arcs);
-	for (size_t i = 0; i < max_num_arcs; ++i)
-		(*arc_prob_cost)[i] = DBL_MAX;
-
-	*arc_fee_cost = tal_arr(ctx, s64, max_num_arcs);
-	for (size_t i = 0; i < max_num_arcs; ++i)
-		(*arc_fee_cost)[i] = INT64_MAX;
-	*arc_capacity = tal_arrz(ctx, s64, max_num_arcs);
-
-	const struct gossmap *gossmap = params->rq->gossmap;
-
-	for (struct gossmap_node *node = gossmap_first_node(gossmap); node;
-	     node = gossmap_next_node(gossmap, node)) {
-		const u32 node_id = gossmap_node_idx(gossmap, node);
-
-		for (size_t j = 0; j < node->num_chans; ++j) {
-			int half;
-			const struct gossmap_chan *c =
-			    gossmap_nth_chan(gossmap, node, j, &half);
-                        struct amount_msat mincap, maxcap;
-
-			if (!channel_is_available(params->rq, c, half))
-				continue;
-
-			/* If a channel cannot forward the total amount we don't
-			 * use it. */
-			if (amount_msat_less(params->amount,
-					     gossmap_chan_htlc_min(c, half)) ||
-			    amount_msat_greater(params->amount,
-						gossmap_chan_htlc_max(c, half)))
-				continue;
-
-			get_constraints(params->rq, c, half, &mincap, &maxcap);
-			/* Assume if min > max, min is wrong */
-			if (amount_msat_greater(mincap, maxcap))
-				mincap = maxcap;
-			/* It is preferable to work on 1msat past the known
-			 * bound. */
-			if (!amount_msat_accumulate(&maxcap, amount_msat(1)))
-				abort();
-
-			/* If amount is greater than the known liquidity upper
-			 * bound we get infinite probability cost. */
-			if (amount_msat_greater_eq(params->amount, maxcap))
-				continue;
-
-			const u32 chan_id = gossmap_chan_idx(gossmap, c);
-
-			const struct gossmap_node *next =
-			    gossmap_nth_node(gossmap, c, !half);
-
-			const u32 next_id = gossmap_node_idx(gossmap, next);
-
-			/* channel to self? */
-			if (node_id == next_id)
-				continue;
-
-			struct arc arc =
-			    arc_from_parts(chan_id, half, 0, false);
-
-			graph_add_arc(*graph, arc, node_obj(node_id),
-				      node_obj(next_id));
-
-			(*arc_capacity)[arc.idx] = 1;
-			(*arc_prob_cost)[arc.idx] =
-			    (-1.0) * log(pickhardt_richter_probability(
-					 mincap, maxcap, params->amount));
-
-			struct amount_msat fee;
-			if (!amount_msat_fee(&fee, params->amount,
-					     c->half[half].base_fee,
-					     c->half[half].proportional_fee))
-				abort();
-			(*arc_fee_cost)[arc.idx] =
-			    fee.millisatoshis + /* Raw: fee co