pylmcf 1.3.0

pylmcf: min-cost-flow solvers for Python and C++

Overview

pylmcf ships two independent deliverables from one source tree:

1. A Python extension — fast min-cost flow via LEMON's solvers, with a stateful object API that warm-restarts on re-solve.
2. A header-only C++ include tree — LEMON extended with warm-restart machinery, plus pylmcf's own solvers: a link-cut tree, two link-cut-tree network simplex variants, and a specialised 1D-chain solver. Downstream C++ code compiles against it directly.

It is used by wnet (efficient Wasserstein and Truncated Wasserstein distance between multidimensional distributions) and wnetalign (alignment of MS or NMR spectra). Both re-solve the same graph hundreds of times with one thing moved, which is what most of the machinery here exists to make cheap.

The LCT simplex variants and 1D chain solver are also available from Python; see the Python solver APIs.

Features

• Fast min-cost flow computation using a C++ backend (LEMON's Network Simplex by default)
• Warm restarts: re-solving after changing supplies, capacities or costs reuses the retained spanning-tree basis instead of rebuilding one — docs
• Multiple solver variants: Network Simplex, Cycle Canceling, Cost Scaling, Capacity Scaling
• Supports capacities, costs, supplies/demands, and per-edge lower bounds (minimum flow)
• NetworkX integration: construct from nx.DiGraph or convert results back for visualization
• Free-threaded CPython support (a separate cp315-abi3t wheel)
• A C++ header tree for downstream packages, including the standalone link-cut-tree data structure

Documentation

See the changelog for release notes and migration details.

Installation

pip install pylmcf

Optional extras for NetworkX support and visualization:

pip install pylmcf[extras]

Usage

Basic usage

import numpy as np
import pylmcf

# 3-node graph with edges 0→1, 0→2, 1→2
# Edges must be sorted by (start, end)
G = pylmcf.Graph(3,
    edge_starts=np.array([0, 0, 1]),
    edge_ends=np.array([1, 2, 2]))

G.set_node_supply(np.array([5, 0, -5]))   # node 0 supplies 5, node 2 demands 5
G.set_edge_costs(np.array([1, 3, 5]))
G.set_edge_capacities(np.array([3, 3, 5]))

G.solve()

G.result()      # np.array([2, 3, 2])  — flow on each edge
G.total_cost()  # 21
G.potentials()  # np.array([-6, -5, 0]) — node potentials (dual solution)

For performance-sensitive dual consumers, use raw_potentials(), dual_values() or allocation-free dual_values_into(). These are also exposed on both LCT variants. See Fast dual certificates.

potentials() returns the dual solution in LEMON's convention: the reduced cost of edge (u, v) is cost + pi[u] - pi[v], and complementary slackness holds against result() (reduced cost > 0 ⇒ flow at its minimum, < 0 ⇒ flow at capacity). Potentials are not unique, so pylmcf returns a canonical choice: of all potentials optimal for the returned flows, the pointwise-largest with pi <= 0 (the smallest with pi >= 0 under "leq" supply, below). Its values are bounded by path costs — LEMON's raw tree potentials can instead carry its internal 2^62 artificial cost on some nodes.

Input arrays may be writable or read-only (including read-only NumPy memory maps); they must be contiguous and aligned to their element size. Returned flow and potential arrays are writable and own their data.

All integer arrays in the OO API are int64, and the dtype is not converted for you — a mismatched array raises TypeError. Costs and minimums must be non-negative, costs at most 2^62 − 1 (network simplex's artificial arcs cost 2^62 and must stay dearer than any real arc), supplies greater than the dtype's minimum (LEMON negates them), and the graph must be feasible (see below for unbalanced supply) or solve() raises RuntimeError: INFEASIBLE. A rejected setter leaves the graph unchanged. As in the functional API below, the value range is unchecked: total supply and every flow must fit in int64 with room to spare, a capacity of INT64_MAX means unbounded, and the total cost must fit in int64.

LCT and 1D-chain solvers

from pylmcf import NetworkSimplexLCT, lmcf_lct, lmcf_lct_dyn, solve_chain_1d

# Same five int64 arrays as the functional LEMON API; arbitrary edge order.
flows = lmcf_lct(supply, starts, ends, capacities, costs)
# lmcf_lct_dyn(...) selects the experimental dynamic-tree variant.
solver = NetworkSimplexLCT(supply, starts, ends, capacities, costs)
solver.solve()
solver.set_node_supply(new_supply)
solver.solve()  # warm restart, with automatic cold fallback

chain = solve_chain_1d(
    np.array([0, 10], dtype=np.int64),
    np.array([5, 0], dtype=np.int64),   # empirical masses
    np.array([0, 5], dtype=np.int64),   # theoretical masses
    kappa=3,
)
chain["total_cost"]  # 15: trash is cheaper than transport

These solvers require contiguous int64 arrays. LCT solvers require balanced supplies, non-negative costs, finite capacities, and zero lower bounds. Costs and topology are fixed on a stateful LCT instance. See Python solver APIs for return values, counters, validation, and the chain's SimpleTrash semantics.

Per-edge lower bounds (minimum flow)

G.set_edge_minimums(np.array([3, 0, 0]))  # edge 0→1 must carry at least 3 units
G.solve()
G.result()      # np.array([3, 2, 3])
G.total_cost()  # 24

Unbalanced supply, and why a problem is infeasible

When total supply and demand differ, the supply type says which side may fall short (LEMON's SupplyType; with balanced supply both mean equality):

G.set_supply_type("geq")  # default: out - in >= supply. Every supply is shipped,
                          # demands may go unmet; needs sum(supply) <= 0.
G.set_supply_type("leq")  # out - in <= supply. Every demand is met, supplies
                          # may go unused; needs sum(supply) >= 0.

Unbalanced re-solves do not warm-restart. When solve() raises INFEASIBLE, infeasibility_cut() says why: a boolean mask of nodes B that must push out more than the edges leaving B can carry (under "leq": must take in more than the edges entering B can deliver). It returns None when the problem is feasible, and does not need or disturb a solve.

G.set_node_supply(np.array([9, 0, -9]))
G.infeasibility_cut()     # array([ True, False, False]) — node 0's out-edges carry at most 6

Re-solving: warm restarts

Changing supplies, capacities or costs and calling solve() again reuses the basis from the previous solve. Nothing to switch on:

H = pylmcf.Graph(3,
    edge_starts=np.array([0, 0, 1]),
    edge_ends=np.array([1, 2, 2]))
H.set_edge_costs(np.array([1, 3, 5]))
H.set_edge_capacities(np.array([3, 3, 5]))

for d in (5, 4, 6, 2):
    H.set_node_supply(np.array([d, 0, -d]))
    H.solve()
    print(d, H.result(), H.total_cost())

H.warm_start_count()    # 2
H.cold_start_count()    # 0
H.dual_repair_count()   # 1

5 [2 3 2] 21
4 [1 3 1] 15
6 [3 3 3] 27
2 [0 2 0] 6

The first solve() is always cold and is not counted. The counters are how you tell whether warm restarts are actually firing — a regression that silently routes every re-solve through a cold rebuild still gives correct answers. Nonzero lower bounds force a cold solve.

The pivot rule (set_pivot_rule(): "block_search" by default, or "first_eligible", "best_eligible", "candidate_list", "altering_list"), the warm-repair strategy (set_warm_repair(): "dual" by default, or "repair_only", "primal", "dual_ratio", "dual_greedy") and the repair time budget (set_warm_repair_budget()) are configurable; the defaults are what you want unless you are benchmarking. Full details, the repair strategies, the set_warm_violation_limit() policy and the C++ API: Warm restarts.

Constructing from a NetworkX graph

Note: Graph.FromNX() iterates over graph elements in Python and is significantly slower than constructing Graph directly from numpy arrays. Prefer the direct API for performance-sensitive code.

import networkx as nx

G_nx = nx.DiGraph()
G_nx.add_edge(0, 1, weight=1, capacity=3)
G_nx.add_edge(0, 2, weight=3, capacity=3)
G_nx.add_edge(1, 2, weight=5, capacity=5)
G_nx.nodes[0]["demand"] = -5
G_nx.nodes[2]["demand"] = 5

G = pylmcf.Graph.FromNX(G_nx)
G.solve()
G.result()      # np.array([2, 3, 2])

Alternative solvers

The default solver is LEMON's Network Simplex. Three alternatives are available via the low-level functional API:

from pylmcf import pylmcf_cpp
import numpy as np

a = lambda x: np.array(x, dtype=np.int64)
supply = a([5, 0, -5])
starts = a([0, 0, 1])
ends   = a([1, 2, 2])
caps   = a([3, 3, 5])
costs  = a([1, 3, 5])

flows = pylmcf_cpp.lmcf(supply, starts, ends, caps, costs)               # Network Simplex (default)
flows = pylmcf_cpp.lmcf_cycle_canceling(supply, starts, ends, caps, costs)
flows = pylmcf_cpp.lmcf_cost_scaling(supply, starts, ends, caps, costs)   # int32/int64 only
flows = pylmcf_cpp.lmcf_capacity_scaling(supply, starts, ends, caps, costs)  # int32/int64 only

Pass return_potentials=True to get (flows, potentials) instead; potentials are always int64 and canonical as above, so every solver returning the same flows returns the same potentials.

The functional API is stateless and has no warm restart. It is duck-typed over int8/16/32/64 (cost and capacity scaling excepted, which need the wider range), and each function has a with- and without-minimums overload.

Value range is the caller's responsibility, at every width. Supplies, capacities, minimums and flows are computed in the input dtype (costs and the objective in int64), and nothing checks that intermediate values fit — that is what keeps the narrow types small. The sum of the positive supplies, the magnitude of the sum of the negative ones, and therefore every arc's flow must fit the dtype with room to spare: at the exact limit some solvers already overflow internally. Outside that range, int8/int16 return wrong results without an error (a feasible problem reported INFEASIBLE, say), and int32/int64 is undefined behaviour. A capacity equal to the dtype's maximum (127 for int8, 32767 for int16, …) means unbounded, not that number. The total cost must fit in int64.

Visualization

G.show()               # display with matplotlib
G.show("graph.png")    # save to file

Threading

Independent Graph objects may be solved concurrently, including on free-threaded CPython — tests/test_free_threading.py hammers 8 threads against a serial oracle. A single Graph must not be driven from two threads at once: solve() mutates the instance, takes no lock, and no binding releases the GIL. See Build modes and threading.

Using the C++ solvers

python -m pylmcf --include

prints the include root; there is nothing to link.

g++ -I$(python -m pylmcf --include) -std=c++20 -O2 mycode.cpp -o mycode

Beyond LEMON, the tree carries pylmcf's own header-only solvers — LinkCutTree, NetworkSimplexLCT, the experimental NetworkSimplexLCTDyn, a lemon::NetworkSimplex-compatible adapter, and ChainSolver1D. None of these have a Python binding. See The C++ header tree.

Requirements

• Python 3.10+
• For a build from source (sdist, or a local checkout): a C++20 compiler

Licence

pylmcf is published under the Boost Software Licence.

src/pylmcf/cpp/lemon is a modified copy of the LEMON graph library, also covered by the Boost Software Licence. lemon/network_simplex.h in particular carries substantial pylmcf-specific work that upstream does not have — the entire warm-restart machinery — so it must not be replaced with a stock LEMON release. Diff, do not overwrite; see the header-tree docs.

Citation

If you use this software, please cite:

Król J, Bochenek M, Jopa S, Kazimierczuk K, Gambin A, Startek MP (2026). WNetAlign: fast and accurate spectra alignment using truncated Wasserstein distance and network simplex. Briefings in Bioinformatics, 27(3), bbag247. https://doi.org/10.1093/bib/bbag247

@article{krol2026wnetalign,
  title   = {WNetAlign: fast and accurate spectra alignment using truncated Wasserstein distance and network simplex},
  author  = {Kr{\'o}l, Justyna and Bochenek, Maria and Jopa, Sylwia and Kazimierczuk, Krzysztof and Gambin, Anna and Startek, Micha{\l} Piotr},
  journal = {Briefings in Bioinformatics},
  volume  = {27},
  number  = {3},
  pages   = {bbag247},
  year    = {2026},
  doi     = {10.1093/bib/bbag247}
}

References

• LEMON Graph Library
• wnet package
• wnetalign package
• Min Cost Flow Problem
• Sleator, D. D. and Tarjan, R. E. (1983). A data structure for dynamic trees. Journal of Computer and System Sciences 26(3), 362–391.
• Cunningham, W. H. (1976). A network simplex method. Mathematical Programming 11(1), 105–116. — the strongly-feasible-basis anti-cycling rule both LCT solvers reproduce.