askrene: add internal API for single-path routes

The single path solver uses the same probability cost and fee cost
estimation of minflow. Single path routes computed this way are
suboptimal with respect to the MCF solution but still are optimal among
any other single path. Computationally is way faster than MCF, therefore
for some trivial payments it should be prefered.

Changelog-None.

Signed-off-by: Lagrang3 <[email protected]>

Author: Lagrang3

plugins/askrene/mcf.c

@@ -2,6 +2,7 @@
 #include <assert.h>
 #include <ccan/asort/asort.h>
 #include <ccan/bitmap/bitmap.h>
+#include <ccan/err/err.h>
 #include <ccan/list/list.h>
 #include <ccan/tal/str/str.h>
 #include <ccan/tal/tal.h>
@@ -159,6 +160,10 @@
  *
  * */
 
+#define PANIC(message)                                                         \
+	errx(1, "Panic in function %s line %d: %s", __func__, __LINE__,        \
+	     message);
+
 #define PARTS_BITS 2
 #define CHANNEL_PARTS (1 << PARTS_BITS)
 
@@ -297,6 +302,10 @@ struct pay_parameters {
 
 	double delay_feefactor;
 	double base_fee_penalty;
+
+	/* keep this here to have a consistent conversion factor all over.
+	 * Probability Cost (x) = prob_cost_factor * (-log Probability(x)) */
+	double prob_cost_factor;
 };
 
 /* Helper function.
@@ -319,6 +328,46 @@ static void set_capacity(s64 *capacity, u64 value, u64 *cap_on_capacity)
 	*cap_on_capacity -= *capacity;
 }
 
+/* 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_sub(&amount, amount, low))
+		return 1.0;
+	if (!amount_msat_sub(&all_states, high, low))
+		PANIC("we expect high > low");
+	if (!amount_msat_sub(&good_states, all_states, amount))
+		PANIC("we expect high > amount");
+	return amount_msat_ratio(good_states, all_states);
+}
+
+/* Returns the probability of success that this channel is able to forward x
+ * considering only the bounds on the liquidity. */
+static double channel_probability(const struct pay_parameters *params,
+				  const struct gossmap_chan *c, const int dir,
+				  struct amount_msat forward)
+{
+	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;
+	/* It is preferable to work on 1msat past the known bound. */
+	if (!amount_msat_accumulate(&maxcap, amount_msat(1)))
+		PANIC("maxcap + 1msat overflows");
+	return pickhardt_richter_probability(mincap, maxcap, forward);
+}
+
 // TODO(eduardo): unit test this
 /* Split a directed channel into parts with linear cost function. */
 static void linearize_channel(const struct pay_parameters *params,
@@ -348,9 +397,8 @@ static void linearize_channel(const struct pay_parameters *params,
 	{
 		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);
+		cost[i] = params->prob_cost_factor * params->cost_fraction[i] /
+			  (b - a);
 	}
 }
 
@@ -867,6 +915,46 @@ get_flow_paths(const tal_t *ctx,
 	return flows;
 }
 
+/* Given a single path build a flow set. */
+static struct flow **
+get_flow_singlepath(const tal_t *ctx, const tal_t *working_ctx,
+		    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 = tal_arr(ctx, struct flow *, 0);
+
+	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;
+	}
+	struct flow *f = tal(ctx, struct flow);
+	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;
+	}
+
+	tal_arr_expand(&flows, f);
+	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
@@ -914,6 +1002,8 @@ struct flow **minflow(const tal_t *ctx,
 			/params->cap_fraction[i];
 	}
 
+	params->prob_cost_factor =
+	    1000 * amount_msat_ratio(params->amount, params->accuracy);
 	params->delay_feefactor = delay_feefactor;
 	params->base_fee_penalty = base_fee_penalty_estimate(amount);
 
@@ -1025,6 +1115,173 @@ static struct amount_msat linear_flows_cost(struct flow **flows,
 	return total;
 }
 
+/* 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);
+
+	/* We use this quantity to compute the cost of routing the entire
+	 * amount, as opposed to the cost per unit of flow. */
+	const u64 pay_amount =
+	    amount_msat_ratio_ceil(params->amount, params->accuracy);
+
+	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);
+
+			if (!gossmap_chan_set(c, half) ||
+			    !c->half[half].enabled)
+				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;
+
+			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] =
+			    params->prob_cost_factor * (-1.0) *
+			    log(channel_probability(params, c, half,
+						    params->amount));
+
+			(*arc_fee_cost)[arc.idx] =
+			    pay_amount *
+			    linear_fee_cost(c->half[half].base_fee,
+					    c->half[half].proportional_fee,
+					    c->half[half].delay,
+					    params->base_fee_penalty,
+					    params->delay_feefactor);
+		}
+	}
+}
+
+/* Similar to minflow but computes routes that have a single path. */
+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_MSAT(1000);
+	params->delay_feefactor = delay_feefactor;
+	params->base_fee_penalty = base_fee_penalty_estimate(amount);
+	params->prob_cost_factor =
+	    1000 * amount_msat_ratio(params->amount, params->accuracy);
+
+	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)) {
+		rq_log(tmpctx, rq, LOG_BROKEN,
+		       "%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, working_ctx, 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;
+}
 
 const char *default_routes(const tal_t *ctx, struct route_query *rq,
 			   const struct gossmap_node *srcnode,

plugins/askrene/mcf.h

@@ -34,6 +34,28 @@ struct flow **minflow(const tal_t *ctx,
 		      double delay_feefactor,
 		      bool single_part);
 
+/**
+ * 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.
+ */
+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);
+
 /* To sanity check: this is the approximation mcf uses for the cost
  * of each channel. */
 struct amount_msat linear_flow_cost(const struct flow *flow,