diff --git a/common/amount.c b/common/amount.c index fe29eef91a71..a8730d2a3d89 100644 --- a/common/amount.c +++ b/common/amount.c @@ -516,6 +516,16 @@ double amount_msat_ratio(struct amount_msat a, struct amount_msat b) return (double)a.millisatoshis / b.millisatoshis; } +u64 amount_msat_ratio_floor(struct amount_msat a, struct amount_msat b) +{ + return a.millisatoshis / b.millisatoshis; +} + +u64 amount_msat_ratio_ceil(struct amount_msat a, struct amount_msat b) +{ + return (a.millisatoshis + b.millisatoshis - 1) / b.millisatoshis; +} + struct amount_msat amount_msat_div(struct amount_msat msat, u64 div) { msat.millisatoshis /= div; diff --git a/common/amount.h b/common/amount.h index 3070038c61f2..dd6ad61bb262 100644 --- a/common/amount.h +++ b/common/amount.h @@ -148,6 +148,12 @@ bool amount_msat_eq_sat(struct amount_msat msat, struct amount_sat sat); /* a / b */ double amount_msat_ratio(struct amount_msat a, struct amount_msat b); +/* floor(a/b) */ +u64 amount_msat_ratio_floor(struct amount_msat a, struct amount_msat b); + +/* ceil(a/b) */ +u64 amount_msat_ratio_ceil(struct amount_msat a, struct amount_msat b); + /* min(a,b) and max(a,b) */ static inline struct amount_msat amount_msat_min(struct amount_msat a, struct amount_msat b) diff --git a/external/Makefile b/external/Makefile index e199108bf0c1..a289032a5f68 100644 --- a/external/Makefile +++ b/external/Makefile @@ -53,6 +53,10 @@ else LDLIBS += -lsodium endif +ifeq ($(HAVE_ZLIB),1) +LDLIBS += -lz +endif + EXTERNAL_LDLIBS := -L${TARGET_DIR} $(patsubst lib%.a,-l%,$(notdir $(EXTERNAL_LIBS))) submodcheck: $(FORCE) diff --git a/plugins/askrene/Makefile b/plugins/askrene/Makefile index 4628e80f2151..a34442df4b27 100644 --- a/plugins/askrene/Makefile +++ b/plugins/askrene/Makefile @@ -1,5 +1,29 @@ -PLUGIN_ASKRENE_SRC := plugins/askrene/askrene.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 -PLUGIN_ASKRENE_HEADER := plugins/askrene/askrene.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 +PLUGIN_ASKRENE_SRC := \ + plugins/askrene/askrene.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_HEADER := \ + plugins/askrene/askrene.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 + PLUGIN_ASKRENE_OBJS := $(PLUGIN_ASKRENE_SRC:.c=.o) $(PLUGIN_ASKRENE_OBJS): $(PLUGIN_ASKRENE_HEADER) @@ -8,3 +32,5 @@ PLUGIN_ALL_SRC += $(PLUGIN_ASKRENE_SRC) PLUGIN_ALL_HEADER += $(PLUGIN_ASKRENE_HEADER) plugins/cln-askrene: $(PLUGIN_ASKRENE_OBJS) $(PLUGIN_LIB_OBJS) $(PLUGIN_COMMON_OBJS) $(JSMN_OBJS) $(CCAN_OBJS) bitcoin/chainparams.o common/gossmap.o common/sciddir_or_pubkey.o common/gossmods_listpeerchannels.o common/fp16.o common/dijkstra.o common/bolt12.o common/bolt12_merkle.o wire/bolt12_wiregen.o wire/onion_wiregen.o common/route.o + +include plugins/askrene/test/Makefile diff --git a/plugins/askrene/algorithm.c b/plugins/askrene/algorithm.c new file mode 100644 index 000000000000..d253325a8b23 --- /dev/null +++ b/plugins/askrene/algorithm.c @@ -0,0 +1,670 @@ +#include "config.h" +#include +#include +#include +#include + +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= 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 new file mode 100644 index 000000000000..40010eb5cfd9 --- /dev/null +++ b/plugins/askrene/algorithm.h @@ -0,0 +1,179 @@ +#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 + +/* 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/graph.c b/plugins/askrene/graph.c new file mode 100644 index 000000000000..973e6adcc365 --- /dev/null +++ b/plugins/askrene/graph.c @@ -0,0 +1,65 @@ +#include "config.h" +#include + +/* 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 new file mode 100644 index 000000000000..e84ed131b938 --- /dev/null +++ b/plugins/askrene/graph.h @@ -0,0 +1,171 @@ +#ifndef LIGHTNING_PLUGINS_ASKRENE_GRAPH_H +#define LIGHTNING_PLUGINS_ASKRENE_GRAPH_H + +/* Defines a graph data structure. */ + +#include "config.h" +#include +#include +#include + +#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< #include #include +#include #include #include #include +#include #include #include #include @@ -163,7 +165,6 @@ static const double CHANNEL_PIVOTS[]={0,0.5,0.8,0.95}; static const s64 INFINITE = INT64_MAX; -static const u32 INVALID_INDEX = 0xffffffff; static const s64 MU_MAX = 100; /* Let's try this encoding of arcs: @@ -234,10 +235,6 @@ static const s64 MU_MAX = 100; * [ 0 1 2 3 4 5 6 ... 31 ] * dual part chandir chanidx */ -struct arc { - u32 idx; -}; - #define ARC_DUAL_BITOFF (0) #define ARC_PART_BITOFF (1) #define ARC_CHANDIR_BITOFF (1 + PARTS_BITS) @@ -289,6 +286,9 @@ struct pay_parameters { // 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]; @@ -304,19 +304,12 @@ struct pay_parameters { * capacity; these quantities remain constant during MCF execution. */ struct linear_network { - u32 *arc_tail_node; - // notice that a tail node is not needed, - // because the tail of arc is the head of dual(arc) - - struct arc *node_adjacency_next_arc; - struct arc *node_adjacency_first_arc; + struct graph *graph; // probability and fee cost associated to an arc double *arc_prob_cost; s64 *arc_fee_cost; s64 *capacity; - - size_t max_num_arcs,max_num_nodes; }; /* This is the structure that keeps track of the network properties while we @@ -330,78 +323,22 @@ struct residual_network { /* potential function on nodes */ s64 *potential; -}; -/* Helper function. - * Given an arc idx, return the dual's idx in the residual network. */ -static struct arc arc_dual(struct arc arc) -{ - arc.idx ^= (1U << ARC_DUAL_BITOFF); - return arc; -} -/* Helper function. */ -static bool arc_is_dual(struct arc arc) -{ - bool dual; - arc_to_parts(arc, NULL, NULL, NULL, &dual); - return dual; -} + /* auxiliary data, the excess of flow on nodes */ + s64 *excess; +}; /* Helper function. * Given an arc of the network (not residual) give me the flow. */ static s64 get_arc_flow( const struct residual_network *network, + const struct graph *graph, const struct arc arc) { - assert(!arc_is_dual(arc)); - assert(arc_dual(arc).idx < tal_count(network->cap)); - return network->cap[ arc_dual(arc).idx ]; -} - -/* Helper function. - * Given an arc idx, return the node from which this arc emanates in the residual network. */ -static u32 arc_tail(const struct linear_network *linear_network, - const struct arc arc) -{ - assert(arc.idx < linear_network->max_num_arcs); - return linear_network->arc_tail_node[ arc.idx ]; -} -/* Helper function. - * Given an arc idx, return the node that this arc is pointing to in the residual network. */ -static u32 arc_head(const struct linear_network *linear_network, - const struct arc arc) -{ - const struct arc dual = arc_dual(arc); - assert(dual.idx < linear_network->max_num_arcs); - return linear_network->arc_tail_node[dual.idx]; -} - -/* Helper function. - * Given node idx `node`, return the idx of the first arc whose tail is `node`. - * */ -static struct arc node_adjacency_begin( - const struct linear_network * linear_network, - const u32 node) -{ - assert(node < linear_network->max_num_nodes); - return linear_network->node_adjacency_first_arc[node]; -} - -/* Helper function. - * Is this the end of the adjacency list. */ -static bool node_adjacency_end(const struct arc arc) -{ - return arc.idx == INVALID_INDEX; -} - -/* Helper function. - * Given node idx `node` and `arc`, returns the idx of the next arc whose tail is `node`. */ -static struct arc node_adjacency_next( - const struct linear_network *linear_network, - const struct arc arc) -{ - assert(arc.idx < linear_network->max_num_arcs); - return linear_network->node_adjacency_next_arc[arc.idx]; + assert(!arc_is_dual(graph, arc)); + struct arc dual = arc_dual(graph, arc); + assert(dual.idx < tal_count(network->cap)); + return network->cap[dual.idx]; } /* Set *capacity to value, up to *cap_on_capacity. Reduce cap_on_capacity */ @@ -426,12 +363,13 @@ static void linearize_channel(const struct pay_parameters *params, if (amount_msat_greater(mincap, maxcap)) mincap = maxcap; - u64 a = mincap.millisatoshis/1000, /* Raw: linearize_channel */ - b = 1 + maxcap.millisatoshis/1000; /* Raw: linearize_channel */ + 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. */ - u64 cap_on_capacity = fp16_to_u64(c->half[dir].htlc_max) / 1000; + u64 cap_on_capacity = + amount_msat_ratio_floor(gossmap_chan_htlc_max(c, dir), params->accuracy); set_capacity(&capacity[0], a, &cap_on_capacity); cost[0]=0; @@ -439,9 +377,9 @@ 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] - *params->amount.millisatoshis /* Raw: linearize_channel */ - /(b-a); + cost[i] = params->cost_fraction[i] * 1000 + * amount_msat_ratio(params->amount, params->accuracy) + / (b - a); } } @@ -456,6 +394,8 @@ alloc_residual_network(const tal_t *ctx, const size_t max_num_nodes, residual_network->cost = tal_arrz(residual_network, s64, max_num_arcs); residual_network->potential = tal_arrz(residual_network, s64, max_num_nodes); + residual_network->excess = + tal_arrz(residual_network, s64, max_num_nodes); return residual_network; } @@ -464,23 +404,25 @@ static void init_residual_network( const struct linear_network * linear_network, struct residual_network* residual_network) { - const size_t max_num_arcs = linear_network->max_num_arcs; - const size_t max_num_nodes = linear_network->max_num_nodes; + const struct graph *graph = linear_network->graph; + const size_t max_num_arcs = graph_max_num_arcs(graph); + const size_t max_num_nodes = graph_max_num_nodes(graph); - for(struct arc arc = {0};arc.idx < max_num_arcs; ++arc.idx) - { - if(arc_is_dual(arc)) + for (struct arc arc = {.idx = 0}; arc.idx < max_num_arcs; ++arc.idx) { + if (arc_is_dual(graph, arc) || !arc_enabled(graph, arc)) continue; - struct arc dual = arc_dual(arc); - residual_network->cap[arc.idx]=linear_network->capacity[arc.idx]; - residual_network->cap[dual.idx]=0; + struct arc dual = arc_dual(graph, arc); + residual_network->cap[arc.idx] = + linear_network->capacity[arc.idx]; + residual_network->cap[dual.idx] = 0; - residual_network->cost[arc.idx]=residual_network->cost[dual.idx]=0; + residual_network->cost[arc.idx] = + residual_network->cost[dual.idx] = 0; } - for(u32 i=0;ipotential[i]=0; + for (u32 i = 0; i < max_num_nodes; ++i) { + residual_network->potential[i] = 0; + residual_network->excess[i] = 0; } } @@ -505,14 +447,17 @@ static int cmp_double(const double *a, const double *b, void *unused) static double get_median_ratio(const tal_t *working_ctx, const struct linear_network* linear_network) { - u64 *u64_arr = tal_arr(working_ctx, u64, linear_network->max_num_arcs/2); - double *double_arr = tal_arr(working_ctx, double, linear_network->max_num_arcs/2); + const struct graph *graph = linear_network->graph; + 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 = {0};arc.idx < linear_network->max_num_arcs; ++arc.idx) { - if (arc_is_dual(arc)) + 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 < linear_network->max_num_arcs/2); + assert(n < max_num_arcs/2); u64_arr[n] = linear_network->arc_fee_cost[arc.idx]; double_arr[n] = linear_network->arc_prob_cost[arc.idx]; n++; @@ -539,10 +484,12 @@ static void combine_cost_function( * Scale by ratio of (positive) medians. */ const double k = get_median_ratio(working_ctx, linear_network); const double ln_30 = log(30); + const struct graph *graph = linear_network->graph; + const size_t max_num_arcs = graph_max_num_arcs(graph); - for(struct arc arc = {0};arc.idx < linear_network->max_num_arcs; ++arc.idx) + for(struct arc arc = {.idx=0};arc.idx < max_num_arcs; ++arc.idx) { - if(arc_tail(linear_network,arc)==INVALID_INDEX) + if (arc_is_dual(graph, arc) || !arc_enabled(graph, arc)) continue; const double pcost = linear_network->arc_prob_cost[arc.idx]; @@ -570,24 +517,12 @@ static void combine_cost_function( } else { residual_network->cost[arc.idx] = combined; } + /* and the respective dual */ + struct arc dual = arc_dual(graph, arc); + residual_network->cost[dual.idx] = -combined; } } -static void linear_network_add_adjacenct_arc( - struct linear_network *linear_network, - const u32 node_idx, - const struct arc arc) -{ - assert(arc.idx < linear_network->max_num_arcs); - linear_network->arc_tail_node[arc.idx] = node_idx; - - assert(node_idx < linear_network->max_num_nodes); - const struct arc first_arc = linear_network->node_adjacency_first_arc[node_idx]; - - linear_network->node_adjacency_next_arc[arc.idx]=first_arc; - linear_network->node_adjacency_first_arc[node_idx]=arc; -} - /* Get the fee cost associated to this directed channel. * Cost is expressed as PPM of the payment. * @@ -597,7 +532,7 @@ static void linear_network_add_adjacenct_arc( * * into * - * fee_microsat = c_fee * x_sat + * fee = c_fee/10^6 * x * * use `base_fee_penalty` to weight the base fee and `delay_feefactor` to * weight the CLTV delay. @@ -626,21 +561,24 @@ struct amount_msat linear_flow_cost(const struct flow *flow, double delay_feefactor) { struct amount_msat msat_cost; - s64 cost = 0; + s64 cost_ppm = 0; double base_fee_penalty = base_fee_penalty_estimate(total_amount); for (size_t i = 0; i < tal_count(flow->path); i++) { const struct half_chan *h = &flow->path[i]->half[flow->dirs[i]]; - cost += linear_fee_cost(h->base_fee, h ->proportional_fee, h->delay, - base_fee_penalty, delay_feefactor); + cost_ppm += + linear_fee_cost(h->base_fee, h->proportional_fee, h->delay, + base_fee_penalty, delay_feefactor); } - - if (!amount_msat_mul(&msat_cost, flow->delivers, cost)) + if (!amount_msat_fee(&msat_cost, flow->delivers, 0, cost_ppm)) abort(); return msat_cost; } +/* FIXME: Instead of mapping one-to-one the indexes in the gossmap, try to + * reduce the number of nodes and arcs used by taking only those that are + * enabled. We might save some cpu if the work with a pruned network. */ static struct linear_network * init_linear_network(const tal_t *ctx, const struct pay_parameters *params) { @@ -651,20 +589,8 @@ init_linear_network(const tal_t *ctx, const struct pay_parameters *params) const size_t max_num_arcs = max_num_chans * ARCS_PER_CHANNEL; const size_t max_num_nodes = gossmap_max_node_idx(gossmap); - linear_network->max_num_arcs = max_num_arcs; - linear_network->max_num_nodes = max_num_nodes; - - linear_network->arc_tail_node = tal_arr(linear_network,u32,max_num_arcs); - for(size_t i=0;iarc_tail_node[i]=INVALID_INDEX; - - linear_network->node_adjacency_next_arc = tal_arr(linear_network,struct arc,max_num_arcs); - for(size_t i=0;inode_adjacency_next_arc[i].idx=INVALID_INDEX; - - linear_network->node_adjacency_first_arc = tal_arr(linear_network,struct arc,max_num_nodes); - for(size_t i=0;inode_adjacency_first_arc[i].idx=INVALID_INDEX; + linear_network->graph = + graph_new(ctx, max_num_nodes, max_num_arcs, ARC_DUAL_BITOFF); linear_network->arc_prob_cost = tal_arr(linear_network,double,max_num_arcs); for(size_t i=0;ihalf[half].base_fee, c->half[half].proportional_fee, @@ -722,25 +649,24 @@ init_linear_network(const tal_t *ctx, const struct pay_parameters *params) // when the `i` hits the `next` node. for(size_t k=0;kgraph, arc, + node_obj(node_id), + node_obj(next_id)); linear_network->capacity[arc.idx] = capacity[k]; linear_network->arc_prob_cost[arc.idx] = prob_cost[k]; - linear_network->arc_fee_cost[arc.idx] = fee_cost; // + the respective dual - struct arc dual = arc_dual(arc); - - linear_network_add_adjacenct_arc(linear_network,next_id,dual); + struct arc dual = arc_dual(linear_network->graph, arc); linear_network->capacity[dual.idx] = 0; linear_network->arc_prob_cost[dual.idx] = -prob_cost[k]; - linear_network->arc_fee_cost[dual.idx] = -fee_cost; } } @@ -749,321 +675,6 @@ init_linear_network(const tal_t *ctx, const struct pay_parameters *params) return linear_network; } -// TODO(eduardo): unit test this -/* Finds an admissible path from source to target, traversing arcs in the - * residual network with capacity greater than 0. - * The path is encoded into prev, which contains the idx of the arcs that are - * traversed. */ - -/* Note we eschew tmpctx here, as this can be called multiple times! */ -static bool -find_admissible_path(const tal_t *working_ctx, - const struct linear_network *linear_network, - const struct residual_network *residual_network, - const u32 source, const u32 target, struct arc *prev) -{ - bool target_found = false; - /* Simple linear queue of node indexes */ - u32 *queue = tal_arr(working_ctx, u32, linear_network->max_num_arcs); - size_t qstart, qend, prev_len = tal_count(prev); - - for(size_t i=0;icap[arc.idx] <= 0) - continue; - - u32 next = arc_head(linear_network,arc); - - assert(next < prev_len); - - // if that node has been seen previously - if(prev[next].idx!=INVALID_INDEX) - continue; - - prev[next] = arc; - assert(qend < linear_network->max_num_arcs); - queue[qend++] = next; - } - } - 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 linear_network* linear_network, - const struct residual_network *residual_network, - const u32 source, - const u32 target, - const struct arc *prev) -{ - s64 flow = INFINITE; - - u32 cur = target; - while(cur!=source) - { - assert(curcap[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(linear_network,arc); - } - - assert(flow0); - return flow; -} - -/* Augment a `flow` amount along the path defined by `prev`.*/ -static void augment_flow( - const struct linear_network *linear_network, - struct residual_network *residual_network, - const u32 source, - const u32 target, - const struct arc *prev, - s64 flow) -{ - u32 cur = target; - - while(cur!=source) - { - assert(cur < tal_count(prev)); - const struct arc arc = prev[cur]; - const struct arc dual = arc_dual(arc); - - assert(arc.idx < tal_count(residual_network->cap)); - assert(dual.idx < tal_count(residual_network->cap)); - - residual_network->cap[arc.idx] -= flow; - residual_network->cap[dual.idx] += flow; - - assert(residual_network->cap[arc.idx] >=0 ); - - // we are traversing in the opposite direction to the flow, - // hence the next node is at the tail of the arc. - cur = arc_tail(linear_network,arc); - } -} - - -// TODO(eduardo): unit test this -/* Finds any flow that satisfy the capacity and balance constraints of the - * uncertainty network. For the balance function condition we have: - * balance(source) = - balance(target) = amount - * balance(node) = 0 , for every other node - * Returns an error code if no feasible flow is found. - * - * 13/04/2023 This implementation uses a simple augmenting path approach. - * */ -static bool find_feasible_flow(const tal_t *working_ctx, - const struct linear_network *linear_network, - struct residual_network *residual_network, - const u32 source, const u32 target, s64 amount) -{ - assert(amount>=0); - - /* path information - * prev: is the id of the arc that lead to the node. */ - struct arc *prev = tal_arr(working_ctx,struct arc,linear_network->max_num_nodes); - - while(amount>0) - { - // find a path from source to target - if (!find_admissible_path(working_ctx, - linear_network, - residual_network, source, target, - prev)) { - return false; - } - - // traverse the path and see how much flow we can send - s64 delta = get_augmenting_flow(linear_network, - residual_network, - source,target,prev); - - // commit that flow to the path - delta = MIN(amount,delta); - assert(delta>0 && delta<=amount); - - augment_flow(linear_network,residual_network,source,target,prev,delta); - amount -= delta; - } - - return true; -} - -// TODO(eduardo): unit test this -/* Similar to `find_admissible_path` but use Dijkstra to optimize the distance - * label. Stops when the target is hit. */ -static bool find_optimal_path(const tal_t *working_ctx, - struct dijkstra *dijkstra, - const struct linear_network *linear_network, - const struct residual_network *residual_network, - const u32 source, const u32 target, - struct arc *prev) -{ - bool target_found = false; - - bitmap *visited = tal_arrz(working_ctx, bitmap, - BITMAP_NWORDS(linear_network->max_num_nodes)); - for(size_t i=0;icap[arc.idx] <= 0) - continue; - - u32 next = arc_head(linear_network,arc); - - s64 cij = residual_network->cost[arc.idx] - - residual_network->potential[cur] - + residual_network->potential[next]; - - // Dijkstra only works with non-negative weights - assert(cij>=0); - - if(distance[next]<=distance[cur]+cij) - continue; - - dijkstra_update(dijkstra,next,distance[cur]+cij); - prev[next]=arc; - } - } - - return target_found; -} - -/* Set zero flow in the residual network. */ -static void zero_flow( - const struct linear_network *linear_network, - struct residual_network *residual_network) -{ - for(u32 node=0;nodemax_num_nodes;++node) - { - residual_network->potential[node]=0; - for(struct arc arc=node_adjacency_begin(linear_network,node); - !node_adjacency_end(arc); - arc = node_adjacency_next(linear_network,arc)) - { - if(arc_is_dual(arc))continue; - - struct arc dual = arc_dual(arc); - - residual_network->cap[arc.idx] = linear_network->capacity[arc.idx]; - residual_network->cap[dual.idx] = 0; - } - } -} - -// TODO(eduardo): unit test this -/* Starting from a feasible flow (satisfies the balance and capacity - * constraints), find a solution that minimizes the network->cost function. - * - * TODO(eduardo) The MCF must be called several times until we get a good - * compromise between fees and probabilities. Instead of re-computing the MCF at - * each step, we might use the previous flow result, which is not optimal in the - * current iteration but I might be not too far from the truth. - * It comes to mind to use cycle cancelling. */ -static bool optimize_mcf(const tal_t *working_ctx, - struct dijkstra *dijkstra, - const struct linear_network *linear_network, - struct residual_network *residual_network, - const u32 source, const u32 target, const s64 amount) -{ - assert(amount>=0); - - zero_flow(linear_network,residual_network); - struct arc *prev = tal_arr(working_ctx,struct arc,linear_network->max_num_nodes); - - const s64 *const distance = dijkstra_distance_data(dijkstra); - - s64 remaining_amount = amount; - - while(remaining_amount>0) - { - if (!find_optimal_path(working_ctx, dijkstra, linear_network, - residual_network, source, target, prev)) { - return false; - } - - // traverse the path and see how much flow we can send - s64 delta = get_augmenting_flow(linear_network,residual_network,source,target,prev); - - // commit that flow to the path - delta = MIN(remaining_amount,delta); - assert(delta>0 && delta<=remaining_amount); - - augment_flow(linear_network,residual_network,source,target,prev,delta); - remaining_amount -= delta; - - // update potentials - for(u32 n=0;nmax_num_nodes;++n) - { - // see page 323 of Ahuja-Magnanti-Orlin - residual_network->potential[n] -= MIN(distance[target],distance[n]); - - /* Notice: - * if node i is permanently labeled we have - * d_i<=d_t - * which implies - * MIN(d_i,d_t) = d_i - * if node i is temporarily labeled we have - * d_i>=d_t - * which implies - * MIN(d_i,d_t) = d_t - * */ - } - } - return true; -} - // flow on directed channels struct chan_flow { @@ -1074,11 +685,11 @@ struct chan_flow * positive balance (returns a node idx with positive balance) * or we discover a cycle (returns a node idx with 0 balance). * */ -static u32 find_path_or_cycle( +static struct node find_path_or_cycle( const tal_t *working_ctx, const struct gossmap *gossmap, const struct chan_flow *chan_flow, - const u32 start_idx, + const struct node source, const s64 *balance, const struct gossmap_chan **prev_chan, @@ -1088,8 +699,8 @@ static u32 find_path_or_cycle( 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 = start_idx; - bitmap_set_bit(visited, start_idx); + 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. */ @@ -1110,9 +721,9 @@ static u32 find_path_or_cycle( /* follow the flow */ if (chan_flow[c_idx].half[dir] > 0) { - const struct gossmap_node *next = + const struct gossmap_node *n = gossmap_nth_node(gossmap, c, !dir); - u32 next_idx = gossmap_node_idx(gossmap, next); + u32 next_idx = gossmap_node_idx(gossmap, n); prev_dir[next_idx] = dir; prev_chan[next_idx] = c; @@ -1135,7 +746,7 @@ static u32 find_path_or_cycle( } bitmap_set_bit(visited, updated_idx); } - return final_idx; + return node_obj(final_idx); } struct list_data @@ -1148,21 +759,23 @@ struct list_data * 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 gossmap *gossmap, - const u32 start_idx, const u32 final_idx, + 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) { - assert(balance[start_idx] < 0); - assert(balance[final_idx] > 0); - s64 delta = -balance[start_idx]; + 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[final_idx]); + 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 = final_idx; cur_idx != start_idx; + 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]; @@ -1183,9 +796,9 @@ static struct flow *substract_flow(const tal_t *ctx, /* Walk again and substract the flow value (delta). */ assert(delta > 0); - balance[start_idx] += delta; - balance[final_idx] -= delta; - for (u32 cur_idx = final_idx; cur_idx != start_idx; + 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]; @@ -1199,12 +812,14 @@ static struct flow *substract_flow(const tal_t *ctx, chan_flow[chan_idx].half[dir] -= delta; } - f->delivers = amount_msat(delta * 1000); + 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 u32 final_idx, +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) @@ -1213,7 +828,7 @@ static void substract_cycle(const struct gossmap *gossmap, const u32 final_idx, u32 cur_idx; /* Compute greatest flow in this cycle. */ - for (cur_idx = final_idx; cur_idx!=INVALID_INDEX;) { + 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); @@ -1221,17 +836,17 @@ static void substract_cycle(const struct gossmap *gossmap, const u32 final_idx, delta = MIN(delta, chan_flow[chan_idx].half[dir]); cur_idx = prev_idx[cur_idx]; - if (cur_idx == final_idx) + if (cur_idx == sink.idx) /* we have come back full circle */ break; } - assert(cur_idx==final_idx); + assert(cur_idx==sink.idx); /* Walk again and substract the flow value (delta). */ assert(delta < INFINITE); assert(delta > 0); - for (cur_idx = final_idx;cur_idx!=INVALID_INDEX;) { + 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); @@ -1239,11 +854,11 @@ static void substract_cycle(const struct gossmap *gossmap, const u32 final_idx, chan_flow[chan_idx].half[dir] -= delta; cur_idx = prev_idx[cur_idx]; - if (cur_idx == final_idx) + if (cur_idx == sink.idx) /* we have come back full circle */ break; } - assert(cur_idx==final_idx); + assert(cur_idx==sink.idx); } /* Given a flow in the residual network, build a set of payment flows in the @@ -1251,16 +866,16 @@ static void substract_cycle(const struct gossmap *gossmap, const u32 final_idx, static struct flow ** get_flow_paths(const tal_t *ctx, const tal_t *working_ctx, - const struct route_query *rq, + const struct pay_parameters *params, const struct linear_network *linear_network, const struct residual_network *residual_network) { struct flow **flows = tal_arr(ctx,struct flow*,0); - const size_t max_num_chans = gossmap_max_chan_idx(rq->gossmap); + 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(rq->gossmap); + 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 @@ -1268,7 +883,7 @@ get_flow_paths(const tal_t *ctx, int *prev_dir = tal_arr(working_ctx,int,max_num_nodes); - u32 *prev_idx = tal_arr(working_ctx,u32,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; @@ -1276,56 +891,52 @@ get_flow_paths(const tal_t *ctx, // Convert the arc based residual network flow into a flow in the // directed channel network. // Compute balance on the nodes. - for(u32 n = 0;ngraph; + 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(linear_network,arc)) + arc = node_adjacency_next(graph,arc)) { - if(arc_is_dual(arc)) + if(arc_is_dual(graph, arc)) continue; - u32 m = arc_head(linear_network,arc); - s64 flow = get_arc_flow(residual_network,arc); + struct node m = arc_head(graph,arc); + s64 flow = get_arc_flow(residual_network, + graph, arc); u32 chanidx; int chandir; - balance[n] -= flow; - balance[m] += flow; + 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(u32 node_idx=0;node_idxgossmap, chan_flow, node_idx, balance, - prev_chan, prev_dir, prev_idx); + while (balance[source.idx] < 0) { + prev_chan[source.idx] = NULL; + struct node sink = find_path_or_cycle( + working_ctx, params->rq->gossmap, chan_flow, source, + balance, prev_chan, prev_dir, prev_idx); - if (balance[final_idx] > 0) + if (balance[sink.idx] > 0) /* case 1. found a path */ { struct flow *fp = substract_flow( - flows, rq->gossmap, node_idx, final_idx, - balance, chan_flow, prev_idx, prev_dir, - prev_chan); + 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(rq->gossmap, final_idx, - chan_flow, prev_idx, prev_dir, - prev_chan); + substract_cycle(params->rq->gossmap, sink, chan_flow, + prev_idx, prev_dir, prev_chan); } } } @@ -1357,12 +968,15 @@ struct flow **minflow(const tal_t *ctx, * 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); - struct dijkstra *dijkstra; params->rq = rq; params->source = source; params->target = target; params->amount = amount; + params->accuracy = AMOUNT_MSAT(1000); + /* FIXME: params->accuracy = amount_msat_max(amount_msat_div(amount, + * 1000), AMOUNT_MSAT(1)); + * */ // template the channel partition into linear arcs params->cap_fraction[0]=0; @@ -1373,9 +987,6 @@ struct flow **minflow(const tal_t *ctx, params->cost_fraction[i]= log((1-CHANNEL_PIVOTS[i-1])/(1-CHANNEL_PIVOTS[i])) /params->cap_fraction[i]; - - // printf("channel part: %ld, fraction: %lf, cost_fraction: %lf\n", - // i,params->cap_fraction[i],params->cost_fraction[i]); } params->delay_feefactor = delay_feefactor; @@ -1383,53 +994,59 @@ struct flow **minflow(const tal_t *ctx, // build the uncertainty network with linearization and residual arcs struct linear_network *linear_network= init_linear_network(working_ctx, params); + const struct graph *graph = linear_network->graph; + const size_t max_num_arcs = graph_max_num_arcs(graph); + const size_t max_num_nodes = graph_max_num_nodes(graph); struct residual_network *residual_network = - alloc_residual_network(working_ctx, linear_network->max_num_nodes, - linear_network->max_num_arcs); - dijkstra = dijkstra_new(working_ctx, gossmap_max_node_idx(rq->gossmap)); + alloc_residual_network(working_ctx, max_num_nodes, max_num_arcs); - const u32 target_idx = gossmap_node_idx(rq->gossmap,target); - const u32 source_idx = gossmap_node_idx(rq->gossmap,source); + const struct node dst = {.idx = gossmap_node_idx(rq->gossmap, target)}; + const struct node src = {.idx = gossmap_node_idx(rq->gossmap, source)}; init_residual_network(linear_network,residual_network); - /* TODO(eduardo): - * Some MCF algorithms' performance depend on the size of maxflow. If we - * were to work in units of msats we 1. risking overflow when computing - * costs and 2. we risk a performance overhead for no good reason. - * - * Working in units of sats was my first choice, but maybe working in - * units of 10, or 100 sats could be even better. - * - * IDEA: define the size of our precision as some parameter got at - * runtime that depends on the size of the payment and adjust the MCF - * accordingly. - * For example if we are trying to pay 1M sats our precision could be - * set to 1000sat, then channels that had capacity for 3M sats become 3k - * flow units. */ - const u64 pay_amount_sats = (params->amount.millisatoshis + 999)/1000; /* Raw: minflow */ - - if (!find_feasible_flow(working_ctx, linear_network, residual_network, - source_idx, target_idx, pay_amount_sats)) { - tal_free(working_ctx); - return NULL; + /* 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, linear_network->graph, src, dst, + residual_network->cap, pay_amount)) { + rq_log(tmpctx, rq, LOG_INFORM, + "%s failed: unable to find a feasible flow.", __func__); + goto fail; } combine_cost_function(working_ctx, linear_network, residual_network, rq->biases, mu); /* We solve a linear MCF problem. */ - if(!optimize_mcf(working_ctx, dijkstra,linear_network,residual_network, - source_idx,target_idx,pay_amount_sats)) - { - tal_free(working_ctx); - return NULL; + if (!mcf_refinement(working_ctx, + linear_network->graph, + residual_network->excess, + residual_network->cap, + residual_network->cost, + residual_network->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, rq, + flow_paths = get_flow_paths(ctx, working_ctx, params, linear_network, residual_network); + 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; } diff --git a/plugins/askrene/priorityqueue.c b/plugins/askrene/priorityqueue.c new file mode 100644 index 000000000000..7f0a96cca699 --- /dev/null +++ b/plugins/askrene/priorityqueue.c @@ -0,0 +1,151 @@ +#define NDEBUG 1 +#include "config.h" +#include + +/* priorityqueue: a data structure for pairs (key, value) with + * 0<=key, 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/priorityqueue.h b/plugins/askrene/priorityqueue.h new file mode 100644 index 000000000000..42d01cc58048 --- /dev/null +++ b/plugins/askrene/priorityqueue.h @@ -0,0 +1,35 @@ +#ifndef LIGHTNING_PLUGINS_ASKRENE_PRIORITYQUEUE_H +#define LIGHTNING_PLUGINS_ASKRENE_PRIORITYQUEUE_H + +/* Defines a priority queue using gheap. */ + +#include "config.h" +#include +#include +#include + +/* 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_PRIORITYQUEUE_H */ diff --git a/plugins/askrene/test/Makefile b/plugins/askrene/test/Makefile new file mode 100644 index 000000000000..88ee0b88a2fd --- /dev/null +++ b/plugins/askrene/test/Makefile @@ -0,0 +1,22 @@ +# Note that these actually #include everything they need, except ccan/ and bitcoin/. +# That allows for unit testing of statics, and special effects. +PLUGIN_ASKRENE_TEST_SRC := $(wildcard plugins/askrene/test/run-*.c) +PLUGIN_ASKRENE_TEST_OBJS := $(PLUGIN_ASKRENE_TEST_SRC:.c=.o) +PLUGIN_ASKRENE_TEST_PROGRAMS := $(PLUGIN_ASKRENE_TEST_OBJS:.o=) + +ALL_C_SOURCES += $(PLUGIN_ASKRENE_TEST_SRC) +ALL_TEST_PROGRAMS += $(PLUGIN_ASKRENE_TEST_PROGRAMS) +$(PLUGIN_RENEPAY_TEST_OBJS): $(PLUGIN_ASKRENE_SRC) + +PLUGIN_ASKRENE_TEST_COMMON_OBJS := + +plugins/askrene/test/run-bfs plugins/askrene/test/run-dijkstra plugins/askrene/test/run-flow plugins/askrene/test/run-mcf plugins/askrene/test/run-mcf-large: \ + plugins/askrene/priorityqueue.o \ + plugins/askrene/graph.o + +$(PLUGIN_ASKRENE_TEST_PROGRAMS): $(PLUGIN_ASKRENE_TEST_COMMON_OBJS) $(PLUGIN_LIB_OBJS) $(PLUGIN_COMMON_OBJS) $(JSMN_OBJS) $(CCAN_OBJS) + +check-askrene: $(PLUGIN_ASKRENE_TEST_PROGRAMS:%=unittest/%) + +check-units: check-askrene + diff --git a/plugins/askrene/test/data/linear_mcf.gz b/plugins/askrene/test/data/linear_mcf.gz new file mode 100644 index 000000000000..b101fd557217 Binary files /dev/null and b/plugins/askrene/test/data/linear_mcf.gz differ diff --git a/plugins/askrene/test/run-bfs.c b/plugins/askrene/test/run-bfs.c new file mode 100644 index 000000000000..5c5f0a045a89 --- /dev/null +++ b/plugins/askrene/test/run-bfs.c @@ -0,0 +1,101 @@ +#include "config.h" +#include +#include +#include +#include +#include +#include + +#define ASKRENE_UNITTEST +#include "../algorithm.c" + +#define MAX_NODES 256 +#define MAX_ARCS 256 +#define DUAL_BIT 7 + +#define CHECK(arg) if(!(arg)){fprintf(stderr, "failed CHECK at line %d: %s\n", __LINE__, #arg); abort();} + +static 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); + } + printf("\n"); +} + +int main(int argc, char *argv[]) +{ + common_setup(argv[0]); + printf("Allocating a memory context\n"); + tal_t *ctx = tal(NULL, tal_t); + assert(ctx); + + printf("Allocating a graph\n"); + struct graph *graph = graph_new(ctx, MAX_NODES, MAX_ARCS, DUAL_BIT); + assert(graph); + + s64 *capacity = tal_arrz(ctx, s64, MAX_ARCS); + struct arc *prev = tal_arr(ctx, struct arc, MAX_NODES); + + graph_add_arc(graph, arc_obj(0), node_obj(1), node_obj(2)); + capacity[0] = 1; + graph_add_arc(graph, arc_obj(1), node_obj(1), node_obj(3)); + capacity[1] = 1; + graph_add_arc(graph, arc_obj(2), node_obj(1), node_obj(6)); + capacity[2] = 1; + graph_add_arc(graph, arc_obj(3), node_obj(2), node_obj(3)); + capacity[3] = 1; + graph_add_arc(graph, arc_obj(4), node_obj(2), node_obj(4)); + capacity[4] = 0; /* disable this arc */ + graph_add_arc(graph, arc_obj(5), node_obj(3), node_obj(4)); + capacity[5] = 1; + graph_add_arc(graph, arc_obj(6), node_obj(3), node_obj(6)); + capacity[6] = 1; + graph_add_arc(graph, arc_obj(7), node_obj(4), node_obj(5)); + capacity[7] = 1; + graph_add_arc(graph, arc_obj(8), node_obj(5), node_obj(6)); + capacity[8] = 1; + + show(graph, node_obj(1)); + show(graph, node_obj(2)); + show(graph, node_obj(3)); + show(graph, node_obj(4)); + show(graph, node_obj(5)); + show(graph, node_obj(6)); + + struct node src = {.idx = 1}; + struct node dst = {.idx = 5}; + + bool result = BFS_path(ctx, graph, src, dst, capacity, 1, prev); + assert(result); + + int pathlen = 0; + int arc_sequence[] = {7, 5, 1}; + int node_sequence[] = {4, 3, 1}; + + printf("path: "); + for (struct node cur = dst; cur.idx != src.idx;) { + struct arc arc = prev[cur.idx]; + printf("node(%" PRIu32 ") arc(%" PRIu32 ") - ", cur.idx, + arc.idx); + cur = arc_tail(graph, arc); + CHECK(pathlen < 3); + CHECK(cur.idx == node_sequence[pathlen]); + CHECK(arc.idx == arc_sequence[pathlen]); + pathlen ++; + } + CHECK(pathlen == 3); + printf("node(%" PRIu32 ") arc(NONE)\n", src.idx); + printf("path length: %d\n", pathlen); + + printf("Freeing memory\n"); + ctx = tal_free(ctx); + + common_shutdown(); + return 0; +} + diff --git a/plugins/askrene/test/run-dijkstra.c b/plugins/askrene/test/run-dijkstra.c new file mode 100644 index 000000000000..3c82e48c03e1 --- /dev/null +++ b/plugins/askrene/test/run-dijkstra.c @@ -0,0 +1,122 @@ +#include "config.h" +#include +#include +#include +#include +#include +#include + +#define ASKRENE_UNITTEST +#include "../algorithm.c" + +// 1->2 7 +// 1->3 9 +// 1->6 14 +// 2->3 10 +// 2->4 15 +// 3->6 2 +// 3->4 11 +// 4->5 6 +// 5->6 9 + +#define MAX_NODES 256 +#define MAX_ARCS 256 +#define DUAL_BIT 7 + +#define CHECK(arg) if(!(arg)){fprintf(stderr, "failed CHECK at line %d: %s\n", __LINE__, #arg); abort();} + +static 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); + } + printf("\n"); +} + +int main(int argc, char *argv[]) +{ + common_setup(argv[0]); + printf("Allocating a memory context\n"); + tal_t *ctx = tal(NULL, tal_t); + assert(ctx); + + printf("Allocating a graph\n"); + struct graph *graph = graph_new(ctx, MAX_NODES, MAX_ARCS, DUAL_BIT); + assert(graph); + + s64 *capacity = tal_arrz(ctx, s64, MAX_ARCS); + s64 *cost = tal_arrz(ctx, s64, MAX_ARCS); + s64 *potential = tal_arrz(ctx, s64, MAX_NODES); + s64 *distance = tal_arr(ctx, s64, MAX_NODES); + struct arc *prev = tal_arr(ctx, struct arc, MAX_NODES); + + graph_add_arc(graph, arc_obj(0), node_obj(1), node_obj(2)); + cost[0] = 7, capacity[0] = 1; + graph_add_arc(graph, arc_obj(1), node_obj(1), node_obj(3)); + cost[1] = 9, capacity[1] = 1; + graph_add_arc(graph, arc_obj(2), node_obj(1), node_obj(6)); + cost[2] = 14, capacity[2] = 1; + graph_add_arc(graph, arc_obj(3), node_obj(2), node_obj(3)); + cost[3] = 10, capacity[3] = 1; + graph_add_arc(graph, arc_obj(4), node_obj(2), node_obj(4)); + cost[4] = 15, capacity[4] = 1; + graph_add_arc(graph, arc_obj(5), node_obj(3), node_obj(4)); + cost[5] = 11, capacity[5] = 1; + graph_add_arc(graph, arc_obj(6), node_obj(3), node_obj(6)); + cost[6] = 2, capacity[6] = 1; + graph_add_arc(graph, arc_obj(7), node_obj(4), node_obj(5)); + cost[7] = 6, capacity[7] = 1; + graph_add_arc(graph, arc_obj(8), node_obj(5), node_obj(6)); + cost[8] = 9, capacity[8] = 1; + + show(graph, node_obj(1)); + show(graph, node_obj(2)); + show(graph, node_obj(3)); + show(graph, node_obj(4)); + show(graph, node_obj(5)); + show(graph, node_obj(6)); + + struct node src = {.idx = 1}; + struct node dst = {.idx = 6}; + + bool result = dijkstra_path(ctx, graph, src, dst, false, capacity, 1, + cost, potential, prev, distance); + CHECK(result); + + int pathlen = 0; + int arc_sequence[] = {6, 1}; + int node_sequence[] = {3, 1}; + + for (struct node cur = dst; cur.idx != src.idx;) { + struct arc arc = prev[cur.idx]; + printf("node(%" PRIu32 ") arc(%" PRIu32 ") - ", cur.idx, + arc.idx); + cur = arc_tail(graph, arc); + CHECK(pathlen < 2); + CHECK(cur.idx == node_sequence[pathlen]); + CHECK(arc.idx == arc_sequence[pathlen]); + pathlen++; + } + CHECK(pathlen == 2); + + for (size_t i = 1; i <= 6; i++) { + printf("node: %zu, distance: %" PRIi64 "\n", i, distance[i]); + } + + CHECK(distance[1] == 0); + CHECK(distance[2] == 7); + CHECK(distance[3] == 9); + CHECK(distance[4] == 20); + CHECK(distance[5] == 26); + CHECK(distance[6] == 11); + + printf("Freeing memory\n"); + ctx = tal_free(ctx); + + common_shutdown(); + return 0; +} diff --git a/plugins/askrene/test/run-flow.c b/plugins/askrene/test/run-flow.c new file mode 100644 index 000000000000..0f1eaba8f8ff --- /dev/null +++ b/plugins/askrene/test/run-flow.c @@ -0,0 +1,114 @@ +#include "config.h" +#include +#include +#include +#include +#include +#include + +#define ASKRENE_UNITTEST +#include "../algorithm.c" + +#define MAX_NODES 256 +#define MAX_ARCS 256 +#define DUAL_BIT 7 + +#define CHECK(arg) if(!(arg)){fprintf(stderr, "failed CHECK at line %d: %s\n", __LINE__, #arg); abort();} + +static void problem1(void){ + printf("Allocating a memory context\n"); + tal_t *ctx = tal(NULL, tal_t); + assert(ctx); + + printf("Allocating a graph\n"); + struct graph *graph = graph_new(ctx, MAX_NODES, MAX_ARCS, DUAL_BIT); + assert(graph); + + s64 *capacity = tal_arrz(ctx, s64, MAX_ARCS); + + graph_add_arc(graph, arc_obj(0), node_obj(1), node_obj(2)); + capacity[0] = 1; + graph_add_arc(graph, arc_obj(1), node_obj(1), node_obj(3)); + capacity[1] = 4; + graph_add_arc(graph, arc_obj(2), node_obj(2), node_obj(4)); + capacity[2] = 1; + graph_add_arc(graph, arc_obj(3), node_obj(2), node_obj(5)); + capacity[3] = 1; + graph_add_arc(graph, arc_obj(4), node_obj(3), node_obj(5)); + capacity[4] = 4; + graph_add_arc(graph, arc_obj(5), node_obj(4), node_obj(6)); + capacity[5] = 1; + graph_add_arc(graph, arc_obj(6), node_obj(6), node_obj(10)); + capacity[6] = 1; + graph_add_arc(graph, arc_obj(7), node_obj(5), node_obj(10)); + capacity[7] = 4; + + struct node src = {.idx = 1}; + struct node dst = {.idx = 10}; + + bool result = simple_feasibleflow(ctx, graph, src, dst, capacity, 5); + CHECK(result); + + CHECK(node_balance(graph, src, capacity) == -5); + CHECK(node_balance(graph, dst, capacity) == 5); + + for (u32 i = 2; i < 10; i++) + CHECK(node_balance(graph, node_obj(i), capacity) == 0); + + printf("Freeing memory\n"); + ctx = tal_free(ctx); +} + +static void problem2(void){ + /* Stress the graph constraints by setting max_num_nodes to exactly the + * number of node that participate and put all nodes in line to achieve + * the largest path length possible. */ + printf("Allocating a memory context\n"); + tal_t *ctx = tal(NULL, tal_t); + assert(ctx); + + printf("Allocating a graph\n"); + struct graph *graph = graph_new(ctx, 5, MAX_ARCS, DUAL_BIT); + assert(graph); + + s64 *capacity = tal_arrz(ctx, s64, MAX_ARCS); + + graph_add_arc(graph, arc_obj(0), node_obj(0), node_obj(1)); + capacity[0] = 1; + graph_add_arc(graph, arc_obj(1), node_obj(1), node_obj(2)); + capacity[1] = 4; + graph_add_arc(graph, arc_obj(2), node_obj(2), node_obj(3)); + capacity[2] = 1; + graph_add_arc(graph, arc_obj(3), node_obj(3), node_obj(4)); + capacity[3] = 1; + + struct node src = {.idx = 0}; + struct node dst = {.idx = 4}; + + bool result = simple_feasibleflow(ctx, graph, src, dst, capacity, 1); + CHECK(result); + + CHECK(node_balance(graph, src, capacity) == -1); + CHECK(node_balance(graph, dst, capacity) == 1); + + for (u32 i = 1; i < 4; i++) + CHECK(node_balance(graph, node_obj(i), capacity) == 0); + + printf("Freeing memory\n"); + ctx = tal_free(ctx); +} + +int main(int argc, char *argv[]) +{ + common_setup(argv[0]); + + printf("\n\nProblem 1\n\n"); + problem1(); + + printf("\n\nProblem 2\n\n"); + problem2(); + + common_shutdown(); + return 0; +} + diff --git a/plugins/askrene/test/run-graph.c b/plugins/askrene/test/run-graph.c new file mode 100644 index 000000000000..c1017e66619e --- /dev/null +++ b/plugins/askrene/test/run-graph.c @@ -0,0 +1,63 @@ +#include "config.h" +#include +#include +#include +#include +#include +#include + +#define ASKRENE_UNITTEST +#include "../graph.c" + +#define MAX_NODES 10 +#define MAX_ARCS 256 +#define DUAL_BIT 7 + +static 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); + } + printf("\n"); +} + +int main(int argc, char *argv[]) +{ + common_setup(argv[0]); + printf("Allocating a memory context\n"); + tal_t *ctx = tal(NULL, tal_t); + assert(ctx); + + printf("Allocating a graph\n"); + struct graph *graph = graph_new(ctx, MAX_NODES, MAX_ARCS, DUAL_BIT); + assert(graph); + + bool success; + + success = graph_add_arc(graph, arc_obj(0), node_obj(0), node_obj(1)); + assert(success); + + success = graph_add_arc(graph, arc_obj(1), node_obj(0), node_obj(2)); + assert(success); + + success = graph_add_arc(graph, arc_obj(2), node_obj(0), node_obj(3)); + assert(success); + + success = graph_add_arc(graph, arc_obj(3), node_obj(3), node_obj(2)); + assert(success); + + show(graph, node_obj(0)); + show(graph, node_obj(1)); + show(graph, node_obj(2)); + show(graph, node_obj(3)); + + printf("Freeing memory\n"); + ctx = tal_free(ctx); + common_shutdown(); + return 0; +} + diff --git a/plugins/askrene/test/run-mcf-large.c b/plugins/askrene/test/run-mcf-large.c new file mode 100644 index 000000000000..ad521956fd59 --- /dev/null +++ b/plugins/askrene/test/run-mcf-large.c @@ -0,0 +1,140 @@ +#include "config.h" +#include +#include +#include +#include +#include +#include + +#define ASKRENE_UNITTEST +#include "../algorithm.c" + +#ifdef HAVE_ZLIB +#include + +gzFile infile; +#define BUFFER_SIZE 10000 +char buffer[BUFFER_SIZE]; + +static int myscanf(const char *fmt, ...) +{ + gzgets(infile, buffer, BUFFER_SIZE); + va_list args; + va_start(args, fmt); + int ret = vsscanf(buffer, fmt, args); + va_end(args); + return ret; +} + +#endif + +#define CHECK(arg) if(!(arg)){fprintf(stderr, "failed CHECK at line %d: %s\n", __LINE__, #arg); abort();} + +/* Read multiple testcases from a file. + * Each test case consist of: + * - one line with N and M, the number of nodes and arcs + * - for each arc a line with numbers head, tail, cap and cost, defining the + * arc's endpoints and values, + * - one line with numbers: amount and best_cost, indicating that we wish to + * send amount from node 0 to node 1 with minimum cost, the correct answer + * should contain a flow cost equal to best_cost. + * + * A feasible solution is guaranteed. + * The last test case has 0 nodes and should be ignored. */ + +static int next_bit(s64 x) +{ + int b; + for (b = 0; (1LL << b) <= x; b++) + ; + return b; +} + +static bool solve_case(const tal_t *ctx) +{ + int ret; + static int c = 0; + c++; + const tal_t *this_ctx = tal(ctx, tal_t); + + int N_nodes, N_arcs; + ret = myscanf("%d %d\n", &N_nodes, &N_arcs); + CHECK(ret == 2); + printf("Testcase %d\n", c); + printf("nodes %d arcs %d\n", N_nodes, N_arcs); + if (N_nodes == 0 && N_arcs == 0) + goto fail; + + const int MAX_NODES = N_nodes; + const int DUAL_BIT = next_bit(N_arcs-1); + const int MAX_ARCS = 1LL << (DUAL_BIT+1); + printf("max nodes %d max arcs %d bit %d\n", MAX_NODES, MAX_ARCS, DUAL_BIT); + + struct graph *graph = graph_new(ctx, MAX_NODES, MAX_ARCS, DUAL_BIT); + CHECK(graph); + + s64 *capacity = tal_arrz(ctx, s64, MAX_ARCS); + s64 *cost = tal_arrz(ctx, s64, MAX_ARCS); + + for (u32 i = 0; i < N_arcs; i++) { + u32 from, to; + ret = myscanf("%" PRIu32 " %" PRIu32 " %" PRIi64 " %" PRIi64, + &from, &to, &capacity[i], &cost[i]); + CHECK(ret == 4); + struct arc arc = {.idx = i}; + graph_add_arc(graph, arc, node_obj(from), node_obj(to)); + + struct arc dual = arc_dual(graph, arc); + cost[dual.idx] = -cost[i]; + } + printf("Reading arcs finished\n"); + struct node src = {.idx = 0}; + struct node dst = {.idx = 1}; + + s64 amount, best_cost; + ret = myscanf("%" PRIi64 " %" PRIi64, &amount, &best_cost); + CHECK(ret == 2); + + bool result = simple_mcf(ctx, graph, src, dst, capacity, amount, cost); + CHECK(result); + + CHECK(node_balance(graph, src, capacity) == -amount); + CHECK(node_balance(graph, dst, capacity) == amount); + + for (u32 i = 2; i < N_nodes; i++) + CHECK(node_balance(graph, node_obj(i), capacity) == 0); + + const s64 total_cost = flow_cost(graph, capacity, cost); + CHECK(total_cost == best_cost); + + tal_free(this_ctx); + return true; + +fail: + tal_free(this_ctx); + return false; +} + +int main(int argc, char *argv[]) +{ +#ifdef HAVE_ZLIB + common_setup(argv[0]); + infile = gzopen("plugins/askrene/test/data/linear_mcf.gz", "r"); + CHECK(infile); + const tal_t *ctx = tal(NULL, tal_t); + CHECK(ctx); + + /* One test case after another. The last test case has N number of nodes + * and arcs equal to 0 and must be ignored. */ + while (solve_case(ctx)) + ; + + ctx = tal_free(ctx); + gzclose(infile); + common_shutdown(); + return 0; +#else + return 0; +#endif +} + diff --git a/plugins/askrene/test/run-mcf.c b/plugins/askrene/test/run-mcf.c new file mode 100644 index 000000000000..aacb51c92ea0 --- /dev/null +++ b/plugins/askrene/test/run-mcf.c @@ -0,0 +1,69 @@ +#include "config.h" +#include +#include +#include +#include +#include +#include + +#define ASKRENE_UNITTEST +#include "../algorithm.c" + +#define CHECK(arg) if(!(arg)){fprintf(stderr, "failed CHECK at line %d: %s\n", __LINE__, #arg); abort();} + +#define MAX_NODES 256 +#define MAX_ARCS 256 +#define DUAL_BIT 7 + +int main(int argc, char *argv[]) +{ + common_setup(argv[0]); + printf("Allocating a memory context\n"); + tal_t *ctx = tal(NULL, tal_t); + assert(ctx); + + printf("Allocating a graph\n"); + struct graph *graph = graph_new(ctx, MAX_NODES, MAX_ARCS, DUAL_BIT); + assert(graph); + + s64 *capacity = tal_arrz(ctx, s64, MAX_ARCS); + s64 *cost = tal_arrz(ctx, s64, MAX_ARCS); + + graph_add_arc(graph, arc_obj(0), node_obj(0), node_obj(1)); + capacity[0] = 2, cost[0] = 0; + graph_add_arc(graph, arc_obj(1), node_obj(0), node_obj(2)); + capacity[1] = 2, cost[1] = 0; + graph_add_arc(graph, arc_obj(2), node_obj(1), node_obj(3)); + capacity[2] = 1, cost[2] = 1; + graph_add_arc(graph, arc_obj(3), node_obj(1), node_obj(4)); + capacity[3] = 1, cost[3] = 2; + graph_add_arc(graph, arc_obj(4), node_obj(2), node_obj(3)); + capacity[4] = 2, cost[4] = 1; + graph_add_arc(graph, arc_obj(5), node_obj(2), node_obj(4)); + capacity[5] = 1, cost[5] = 2; + graph_add_arc(graph, arc_obj(6), node_obj(3), node_obj(5)); + capacity[6] = 3, cost[6] = 0; + graph_add_arc(graph, arc_obj(7), node_obj(4), node_obj(5)); + capacity[7] = 3, cost[7] = 0; + + struct node src = {.idx = 0}; + struct node dst = {.idx = 5}; + + bool result = simple_mcf(ctx, graph, src, dst, capacity, 4, cost); + CHECK(result); + + CHECK(node_balance(graph, src, capacity) == -4); + CHECK(node_balance(graph, dst, capacity) == 4); + + for (u32 i = 1; i < 4; i++) + CHECK(node_balance(graph, node_obj(i), capacity) == 0); + + const s64 total_cost = flow_cost(graph, capacity, cost); + printf("best flow cost: %" PRIi64 "\n", total_cost); + CHECK(total_cost == 5); + + printf("Freeing memory\n"); + ctx = tal_free(ctx); + common_shutdown(); + return 0; +} diff --git a/plugins/askrene/test/run-pqueue.c b/plugins/askrene/test/run-pqueue.c new file mode 100644 index 000000000000..bd7850b935a3 --- /dev/null +++ b/plugins/askrene/test/run-pqueue.c @@ -0,0 +1,105 @@ +#include "config.h" +#include +#include +#include +#include +#include + +#define ASKRENE_UNITTEST +#include "../priorityqueue.c" + +#define CHECK(arg) if(!(arg)){fprintf(stderr, "failed CHECK at line %d: %s\n", __LINE__, #arg); abort();} + +static void priorityqueue_show(struct priorityqueue *q) +{ + printf("size of queue: %zu\n", priorityqueue_size(q)); + printf("empty?: %s\n", priorityqueue_empty(q) ? "true" : "false"); + if (!priorityqueue_empty(q)) + printf("top of the queue: %" PRIu32 "\n", priorityqueue_top(q)); + const s64 *value = priorityqueue_value(q); + for (u32 i = 0; i < priorityqueue_maxsize(q); i++) { + printf("(%" PRIu32 ", %" PRIi64 ")", i, value[i]); + } + + printf("\n\n"); +} + +int main(int argc, char *argv[]) +{ + common_setup(argv[0]); + printf("Hello world!\n"); + + printf("Allocating a memory context\n"); + tal_t *ctx = tal(NULL, tal_t); + CHECK(ctx); + + printf("Allocating a priorityqueue\n"); + struct priorityqueue *q; + q = priorityqueue_new(ctx, 5); + CHECK(q); + + /* reset all values */ + priorityqueue_init(q); + priorityqueue_show(q); + CHECK(priorityqueue_empty(q)); + CHECK(priorityqueue_size(q)==0); + + priorityqueue_update(q, 0, 10); + priorityqueue_show(q); + CHECK(!priorityqueue_empty(q)); + CHECK(priorityqueue_size(q)==1); + CHECK(priorityqueue_top(q)==0); + + priorityqueue_update(q, 0, 3); + priorityqueue_show(q); + CHECK(!priorityqueue_empty(q)); + CHECK(priorityqueue_size(q)==1); + CHECK(priorityqueue_top(q)==0); + + priorityqueue_update(q, 1, 3); + priorityqueue_show(q); + CHECK(!priorityqueue_empty(q)); + CHECK(priorityqueue_size(q)==2); + // CHECK(priorityqueue_top(q)==0); + + priorityqueue_update(q, 1, 5); + priorityqueue_show(q); + CHECK(!priorityqueue_empty(q)); + CHECK(priorityqueue_size(q)==2); + CHECK(priorityqueue_top(q)==0); + + priorityqueue_update(q, 1, -1); + priorityqueue_show(q); + CHECK(!priorityqueue_empty(q)); + CHECK(priorityqueue_size(q)==2); + CHECK(priorityqueue_top(q)==1); + + priorityqueue_pop(q); + priorityqueue_show(q); + CHECK(!priorityqueue_empty(q)); + CHECK(priorityqueue_size(q)==1); + CHECK(priorityqueue_top(q)==0); + + priorityqueue_update(q, 1, 0); + priorityqueue_show(q); + CHECK(!priorityqueue_empty(q)); + CHECK(priorityqueue_size(q)==2); + CHECK(priorityqueue_top(q)==1); + + priorityqueue_update(q, 4, -10); + priorityqueue_show(q); + CHECK(!priorityqueue_empty(q)); + CHECK(priorityqueue_size(q)==3); + CHECK(priorityqueue_top(q)==4); + + priorityqueue_pop(q); + priorityqueue_show(q); + CHECK(!priorityqueue_empty(q)); + CHECK(priorityqueue_size(q)==2); + CHECK(priorityqueue_top(q)==1); + + printf("Freeing memory\n"); + ctx = tal_free(ctx); + common_shutdown(); + return 0; +} diff --git a/tests/test_askrene.py b/tests/test_askrene.py index 2607beb6f384..580215a70653 100644 --- a/tests/test_askrene.py +++ b/tests/test_askrene.py @@ -557,16 +557,21 @@ def test_getroutes_fee_fallback(node_factory): # 0 -> 1 -> 3: high capacity, high fee (1%) # 0 -> 2 -> 3: low capacity, low fee. + # (We disable reverse, since it breaks median calc!) gsfile, nodemap = generate_gossip_store([GenChannel(0, 1, capacity_sats=20000, - forward=GenChannel.Half(propfee=10000)), + forward=GenChannel.Half(propfee=10000), + reverse=GenChannel.Half(enabled=False)), GenChannel(0, 2, - capacity_sats=10000), + capacity_sats=10000, + reverse=GenChannel.Half(enabled=False)), GenChannel(1, 3, capacity_sats=20000, - forward=GenChannel.Half(propfee=10000)), + forward=GenChannel.Half(propfee=10000), + reverse=GenChannel.Half(enabled=False)), GenChannel(2, 3, - capacity_sats=10000)]) + capacity_sats=10000, + reverse=GenChannel.Half(enabled=False))]) # Set up l1 with this as the gossip_store l1 = node_factory.get_node(gossip_store_file=gsfile.name) @@ -1037,10 +1042,10 @@ def test_real_data(node_factory, bitcoind): # CI, it's slow. if SLOW_MACHINE: limit = 25 - expected = (4, 25, 1533317, 143026, 91) + expected = (6, 25, 1544756, 142986, 91) else: limit = 100 - expected = (8, 95, 6007785, 564997, 91) + expected = (9, 95, 6347877, 566288, 92) fees = {} for n in range(0, limit): @@ -1154,10 +1159,10 @@ def test_real_biases(node_factory, bitcoind): # CI, it's slow. if SLOW_MACHINE: limit = 25 - expected = ({1: 4, 2: 5, 4: 7, 8: 11, 16: 14, 32: 19, 64: 25, 100: 25}, 0) + expected = ({1: 5, 2: 7, 4: 7, 8: 11, 16: 14, 32: 19, 64: 25, 100: 25}, 0) else: limit = 100 - expected = ({1: 19, 2: 25, 4: 36, 8: 51, 16: 66, 32: 81, 64: 96, 100: 96}, 0) + expected = ({1: 23, 2: 31, 4: 40, 8: 53, 16: 70, 32: 82, 64: 96, 100: 96}, 0) l1.rpc.askrene_create_layer('biases') num_changed = {} @@ -1201,8 +1206,15 @@ def amount_through_chan(chan, routes): if route2 != route: # It should have avoided biassed channel amount_after = amount_through_chan(chan, route2['routes']) - assert amount_after < amount_before - num_changed[bias] += 1 + if amount_after < amount_before: + num_changed[bias] += 1 + else: + # We bias -4 against 83x88x31908/0 going to node 83, and this is violated. + # Both routes contain three paths, all via 83x88x31908/0. + # The first amounts 49490584, 1018832, 49490584, + # The second amounts 25254708, 25254708, 49490584, + # Due to fees and rounding, we actually spend 1msat more on the second case! + assert (n, bias, chan) == (83, 4, '83x88x31908/0') # Undo bias l1.rpc.askrene_bias_channel(layer='biases', short_channel_id_dir=chan, bias=0)