mixle.relations module

Relations over structured spaces, enumerated in order of a residual.

A Relation is not an optimization problem – it is a constraint imposed on a structured space (matchings, spanning trees, strings near a center, hidden-state sequences, feature subsets), whose members are enumerated ranked by a residual/cost. Finding the single best member is incidental; the value is the whole ranked set. You specify the relation, then ask it for an enumerator() – the same shape as a distribution yielding a sampler() / estimator() / enumerator(). Every relation shares one surface:

relation.solve()       -> the minimal-residual Solution (or None if the relation is empty)
relation.top(k)        -> the k smallest-residual members as a list
relation.enumerator()  -> a lazy iterator over members, smallest residual first
for solution in relation: ...

Each item is a Solution namedtuple (value, objective) – it reads as sol.value / sol.objective and still unpacks as value, objective = sol. value is the member itself (an assignment, a nearby string, a state sequence, a feature subset, …) and objective is its residual: a cost (minimized) or score (maximized); sense records which.

Assignment, spanning tree, the edit-distance ball, k-best Viterbi, shortest path, and best-subset regression are all specified and consumed the same way, each delegating to whatever engine fits (Murty for assignment, Gabow for spanning trees, A* / best_first_paths() for paths, mixle.enumeration.hmm_paths.hmm_best_paths() for Viterbi, Dijkstra / nearest_first() for the edit-distance ball, exhaustive ranking for best-subset). The two shared low-level engines are best_first_paths() (k-best paths to a goal) and nearest_first() (distinct states outward from a center – an expanding metric ball).

>>> from mixle.relations import Assignment
>>> sol = Assignment([[1, 9], [9, 1]]).solve()
>>> sol.value, sol.objective          # the column assignment and its total cost
(array([0, 1]), 2.0)
admm_bounded_least_squares(a, b, lower=0.0, upper=np.inf, *, rho=1.0, max_iter=5000, tol=1.0e-8)[source]

Solve min_x ||A x - b||^2 subject to lower <= x <= upper by ADMM.

The alternating-direction method of multipliers splits the problem as f(x) = ||A x - b||^2 plus the box indicator g(z), with x = z, and alternates: an x-update (the ridge solve (A^T A + rho I) x = A^T b + rho (z - u), factorized once), a z-update (project x + u onto the box), and the scaled dual update u += x - z. This is the augmented-Lagrangian path “beyond pure penalty”: it converges to the exact constrained optimum (lower=0, upper=inf recovers non-negative least squares). Returns the bounded solution x.

lower <= upper elementwise, neither NaN; rho > 0 (the ridge term is SPD only then); max_iter >= 1. An infinite upper stays legal – it is how the NNLS case is spelled.

Parameters:
Return type:

ndarray

class Assignment(cost, maximize=False)[source]

Bases: Relation

Linear assignment / bipartite matching: match rows to columns at extremal total cost.

The solution value is col_ind – col_ind[i] is the column assigned to row i.

Parameters:

  • cost (np.ndarray)

  • maximize (bool)

enumerator(k=None)[source]

Yield up to k assignments ordered by total assignment cost.

Parameters:

k (int | None)

Return type:

Iterator[Solution]

class BestSubsetRegression(X, y, criterion='aic', max_size=None, intercept=True)[source]

Bases: Relation

Best-subset feature selection for least squares, enumerated in increasing selection criterion.

Solution values are feature-index tuples ranked by criterion: residual sum of squares ("rss"), Akaike ("aic") or Bayesian ("bic") information criterion (Gaussian form). Best-subset selection is inherently exponential, so this scores subsets exhaustively up to max_size features – cap max_size (and/or the number of features) for large p.

Parameters:

  • X (np.ndarray) – Design matrix, shape (n, p).

  • y (np.ndarray) – Response vector, length n.

  • criterion (str) – "aic" (default), "bic", or "rss".

  • max_size (int | None) – Largest subset size to consider (None = all p features).

  • intercept (bool) – Fit an (unpenalized, always-included) intercept column.

enumerator(k=None)[source]

Yield up to k feature subsets ordered by the configured criterion.

Parameters:

k (int | None)

Return type:

Iterator[Solution]

branch_and_bound_milp(c, a_ub=None, b_ub=None, integer=None, bounds=None, *, sense='min', tol=1.0e-6)[source]

Solve a mixed-integer linear program by branch-and-bound over the LP relaxation.

Minimizes (sense="min") or maximizes (sense="max") c @ x subject to a_ub @ x <= b_ub and per-variable bounds (lo, hi), with the variables indexed by integer constrained to integers (default: all). Returns (objective, x) or None if infeasible. Each node solves the continuous relaxation with scipy.optimize.linprog (HiGHS) and, if an integer variable is fractional, branches into x_i <= floor and x_i >= ceil; best-bound search prunes nodes that cannot beat the incumbent. Exact for bounded integer feasible regions.

tol is the integrality/pruning tolerance and must be finite and strictly below 0.5: it is the width within which a relaxation coordinate counts as already integral, and half a unit is the most that can ever mean. See the note at its validation below for what a wider one returned.

Parameters:
Return type:

tuple[float, ndarray] | None

cardinality_constrained_milp(c, a_ub, b_ub, max_nonzero, bounds, *, sense='min')[source]

Minimize/maximize c @ x with at most max_nonzero of the variables nonzero.

Adds a cardinality (sparsity) constraint to the linear program a_ub @ x <= b_ub with per-variable bounds via the standard big-M indicator formulation: a binary z_i gates each variable (lower_i z_i <= x_i <= upper_i z_i, so z_i = 0 forces x_i = 0) and sum z_i <= max_nonzero; the extended mixed-integer program is solved by branch_and_bound_milp(). Every bound must be finite – the bounds are the big-M constants of the indicator rows. Returns (objective, x) (the sparse optimizer) or None if infeasible. This is the indicator/ set-membership/cardinality constraint primitive (best-subset selection, sparse design).

Parameters:
Return type:

tuple[float, ndarray] | None

class Design(cost, open, flow)[source]

Bases: NamedTuple

A network-design solution: cost (total), open (opened-arc/facility mask), and the induced flow.

Parameters:
cost: float

Alias for field number 0

open: ndarray

Alias for field number 1

flow: ndarray

Alias for field number 2

class EditDistance(center, alphabet, *, sub_cost=None, ins_cost=None, del_cost=None, max_distance=None)[source]

Bases: Relation

Enumerate strings outward from a center by (non-uniform) edit distance – an edit-distance ball.

You give a single center string and an alphabet, not two endpoints (the distance between two fixed strings is just one number). The enumerator yields strings in increasing edit distance from the center: the center itself at distance 0, then its 1-edit neighbours, then 2-edit, and so on – a Dijkstra expansion over string space with per-operation costs (nearest_first()). The ball is infinite (insertions grow strings without bound), so bound it with max_distance or by taking top(k) / enumerator(k). Solution values are strings (or symbol tuples, matching the center’s type); the objective is the edit distance from the center.

Parameters:

  • center (Iterable[Any]) – The center string (or sequence of symbols).

  • alphabet (Iterable[Any]) – The symbols available for substitution and insertion.

  • sub_cost (Callable[[Any, Any], float] | None) – (a, b) -> cost of substituting a with b (default unit; 0 if equal).

  • ins_cost (Callable[[Any], float] | None) – c -> cost of inserting symbol c (default 1).

  • del_cost (Callable[[Any], float] | None) – a -> cost of deleting symbol a (default 1).

  • max_distance (float | None) – Only enumerate strings within this edit distance (None = unbounded/lazy).

enumerator(k=None)[source]

Yield edit-distance neighbors ordered outward from the center.

Parameters:

k (int | None)

Return type:

Iterator[Solution]

class Flow(value, flow)[source]

Bases: NamedTuple

A resolved flow: value (total/objective cost) and flow[u, v] per-arc, like max_flow’s output.

Parameters:
value: float

Alias for field number 0

flow: ndarray

Alias for field number 1

graph_coloring(adjacency)[source]

Exact minimum proper vertex coloring of an undirected graph.

adjacency is an n x n symmetric 0/1 (or boolean) matrix with a zero diagonal. Returns (k, coloring) where k is the chromatic number and coloring[v] in 0..k-1 gives no two adjacent vertices the same color. Solved by backtracking with the standard symmetry break (a vertex may introduce at most one new color), trying k = 1, 2, ... until colorable – exact, but worst-case exponential, so intended for small/medium graphs.

Parameters:

adjacency (Any)

Return type:

tuple[int, list[int]]

irreducible_infeasible_subset(a_ub, b_ub, bounds=None)[source]

Find an irreducible infeasible subset (IIS) of the linear constraints a_ub @ x <= b_ub.

Returns the row indices of a minimal infeasible subset: the subsystem is itself infeasible, yet dropping any single one of its rows makes it feasible (within the variable bounds, default unbounded). Returns None if the full system is already feasible. Uses the deletion filter – tentatively remove each constraint and keep it removed whenever the remainder stays infeasible – so the result certifies which constraints conflict, the standard infeasibility diagnostic.

Parameters:
Return type:

list[int] | None

min_cost_flow(cap, cost, supply)[source]

Min-cost feasible flow meeting node supply (positive = source) under arc cap/cost.

Successive-shortest-path with node potentials: a super-source feeds every supply > 0 node and every supply < 0 node feeds a super-sink (arc capacity |supply|, cost 0), then the residual super-source -> super-sink path of least reduced cost is repeatedly augmented (Dijkstra, since the potentials keep every reduced cost non-negative) until the super-sink is unreachable. Potentials are seeded with a Bellman-Ford pass so arbitrary-sign cost entries are handled, then updated by the per-round Dijkstra distances (the standard Johnson’s-algorithm maintenance).

Successive-shortest-path alone only guarantees minimality among source-to-sink augmenting paths – it says nothing about a profitable cycle that never touches the super-source/super-sink at all. Once feasibility is reached, _cancel_negative_cycles() runs a Bellman-Ford negative-cycle-canceling pass over the resulting residual graph and pushes flow around any negative-cost cycle it finds (up to its bottleneck capacity) – which strictly lowers cost without disturbing any node’s supply/demand balance, since a cycle has zero net flow at every node it passes through. This is what makes the returned flow truly minimum rather than merely feasible.

cap/cost are n x n arc matrices; supply is length n and must sum to (approximately) zero. Returns Flow(value=total cost, flow=(n, n) arc flows); raises ValueError if the supply cannot be fully routed under the given capacities, or if the cost structure admits an unbounded negative cycle (infinite capacity available to an improving cycle, so no finite minimum cost exists).

Parameters:
Return type:

Flow

multicommodity_flow(cap, cost, demands)[source]

Multi-commodity min-cost flow (different products/grades share arc cap) via LP relaxation.

cap/cost are n x n arc matrices (an arc exists where cap > 0); demands lists (source, sink, amount) rows, one per commodity. Each commodity gets its own flow variables with its own conservation constraints (per-commodity supply at its source, demand at its sink, balance elsewhere); every arc’s combined commodity flow is capped at cap[u, v]. Solved as a single LP (scipy.optimize.linprog, HiGHS) over the stacked per-commodity arc variables. Returns Flow with the aggregate cost and the summed arc flows (flow = sum_k x_k); raises ValueError if the demands cannot be met under the shared capacities.

Parameters:
Return type:

Flow

network_design(nodes, arcs, fixed_costs, demands)[source]

Capacitated fixed-charge network design (which arcs/facilities to open) via big-M MILP.

nodes lists node ids; arcs lists candidate directed arcs (fixed_costs[i] is the cost of opening arcs[i]); demands is a length-len(nodes) net supply/demand vector aligned with nodes (positive = supply, negative = demand, summing to zero) – the same convention as min_cost_flow()’s supply. Continuous flow x[arc] >= 0 is big-M-linked to the binary open/closed decision y[arc] (x[arc] <= M * y[arc], M sized from the total demand so it never binds a genuinely open arc); flow balance A x = demands is encoded as the paired inequalities A x <= demands and -A x <= -demands (branch_and_bound_milp() exposes only a_ub/b_ub). The frozen signature carries no per-unit routing cost, so the objective is the total opening cost fixed_costs @ y (the classic uncapacitated-routing fixed-charge network design: which arcs to build so the demand vector is routable at minimum build cost). Returns Design(cost, open, flow) where open is the boolean opened-arc mask and flow is the induced (n, n) arc flow; raises ValueError if no arc-opening choice can route the demands.

Parameters:
Return type:

Design

class Relation[source]

Bases: ABC

A constraint over a structured space whose members are enumerated ranked by a residual.

Subclasses implement enumerator() (yielding Solution items); solve(), top() and iteration come for free. sense is "min" (residual minimized, members out in increasing cost) or "max" (residual maximized, members out in decreasing score).

abstractmethod enumerator(k=None)[source]

Lazily yield Solution items best-first; at most k if given (None = all).

Parameters:

k (int | None)

Return type:

Iterator[Solution]

solve()[source]

The single optimal Solution, or None if the problem is infeasible.

Return type:

Solution | None

top(k)[source]

The k best solutions as a list.

k is an exact non-negative count: top(0) used to return ONE solution and top(2.9) three, because the engines counted emissions after yielding and compared against the raw float (MXR-080-1903).

Parameters:

k (int)

Return type:

list[Solution]

sampler(seed=None, *, temperature=1.0, k=None, uniform=False, rng=None)[source]

Return a RelationSampler that draws members under a Gibbs measure over the objective.

A relation is a specification of a structured space, not itself a random object; sampling it means imposing a distribution over its members, which needs an RNG and a temperature. So – like every other mixle object – it hands back a sampler (relation.sampler(seed).sample(size)) that owns the stream and the Gibbs measure, rather than being sampled directly.

Each enumerated member is weighted exp(-objective / temperature) when sense == "min" (low cost favoured) or exp(objective / temperature) when sense == "max". temperature -> 0 concentrates on the optimum; -> inf (or uniform=True) is uniform. The draw is an exact Gibbs sample only when the relation is finite and fully enumerated (k=None); pass k to truncate an infinite/large relation to its k best (the dropped tail is the lowest-weight mass – a good low-temperature approximation, and k is required if infinite).

Parameters:

  • seed (int | None) – scalar seed for the sampler’s RandomState (ignored if rng is given).

  • temperature (float) – Gibbs temperature (default 1.0).

  • k (int | None) – enumerate at most this many members (None = all; required if infinite).

  • uniform (bool) – ignore objectives and sample uniformly over the enumerated members.

  • rng – a shared numpy.random.RandomState (takes precedence over seed).

Return type:

RelationSampler

class RelationSampler(relation, seed=None, *, temperature=1.0, k=None, uniform=False, rng=None)[source]

Bases: object

Draws members of a Relation under a Gibbs measure exp(-objective / temperature).

Constructed via Relation.sampler(). It enumerates the relation’s members once (lazily, on the first draw) and caches the resulting categorical, so repeated sample calls are low-overhead. size=None returns one member value; size=int returns a list of that many draws.

Parameters:

  • relation (Relation)

  • seed (int | None)

  • temperature (float)

  • k (int | None)

  • uniform (bool)

sample(size=None)[source]

Draw a member value (size=None) or a list of size member values.

Parameters:

size (int | None)

Return type:

Any

class ShortestPath(start, successors, is_goal=None, *, sense='min', heuristic=None)[source]

Bases: Relation

k-shortest-path / best-first search over an arbitrary state graph.

Specify the graph by start and successors; the solution value is the list of states from start to a goal. By default a sink (a state with no successors) is a goal, so a finite DAG search needs no goal test; pass is_goal for infinite graphs or early goals. Use sense="max" for highest-score paths and an admissible heuristic for A*.

Parameters:

  • start (Any)

  • successors (Callable[[Any], Iterable[tuple[Any, float]]])

  • is_goal (Callable[[Any], bool] | None)

  • sense (str)

  • heuristic (Callable[[Any], float] | None)

enumerator(k=None)[source]

Yield up to k shortest or highest-scoring paths.

Parameters:

k (int | None)

Return type:

Iterator[Solution]

max_clique(adjacency)[source]

A maximum clique (largest mutually-adjacent vertex set) of an undirected graph.

adjacency is an n x n symmetric 0/1 (or boolean) matrix with a zero diagonal. Returns the sorted vertices of one maximum clique via Carraghan-Pardalos branch-and-bound (prune when the current clique plus the remaining candidates cannot beat the incumbent) – exact, worst-case exponential, intended for small/medium graphs.

Parameters:

adjacency (Any)

Return type:

list[int]

max_independent_set(adjacency)[source]

A maximum independent set (largest pairwise-non-adjacent vertex set) – a max clique of the complement.

Parameters:

adjacency (Any)

Return type:

list[int]

is_stable_matching(match, proposer_prefs, receiver_prefs)[source]

Return True iff match has no blocking pair (a mutually-preferred unmatched proposer/receiver).

Parameters:
Return type:

bool

max_flow(capacity, source, sink)[source]

Maximum source -> sink flow in a directed network (Edmonds-Karp).

capacity is an n x n non-negative matrix of arc capacities. Returns (value, flow) where flow[u, v] is the flow on arc u -> v (conserved at every node but the source/sink) and value is the total flow out of source. Edmonds-Karp augments along BFS shortest paths in the residual network, so it runs in O(V E^2) and terminates on real-valued capacities.

capacity must be square and NaN-free, and source/sink must be distinct in-range nodes; see _require_square(), _reject_nan() and the note below for what each of those used to do instead (MXR-080-1903).

Parameters:
Return type:

tuple[float, ndarray]

min_arborescence(weight, root=0)[source]

Minimum-weight spanning arborescence rooted at root (directed MST; Chu-Liu/Edmonds).

weight is an n x n matrix of directed arc costs with +inf for absent arcs. Returns (total, parent) where parent[v] is the chosen in-arc tail for each non-root v (and parent[root] = -1), forming the minimum-cost arborescence in which every node is reachable from root; returns None if no such arborescence exists.

NaN and -inf are refused for the same reason as in tsp_held_karp(): arc presence is decided by np.isfinite, so both used to read as “absent arc” (MXR-080-1903).

Parameters:
Return type:

tuple[float, list[int]] | None

min_cut(capacity, source, sink)[source]

Minimum source/sink cut of a directed network (via max-flow; the max-flow min-cut theorem).

Returns (capacity, source_side, cut_edges): the cut capacity (equal to the max-flow value), the set of nodes on the source side (reachable from source in the final residual graph), and the saturated arcs crossing from the source side to the sink side.

Input contract as max_flow()’s, and enforced there: a NaN capacity used to make this return capacity 0.0 with an EMPTY cut-edge list on a graph that plainly has arcs (MXR-080-1903).

Parameters:
Return type:

tuple[float, list[int], list[tuple[int, int]]]

stable_matching(proposer_prefs, receiver_prefs)[source]

Proposer-optimal stable matching via Gale-Shapley.

proposer_prefs[i] is proposer i’s receivers in descending preference; receiver_prefs[j] likewise for receiver j. Preference lists may be partial (an unlisted partner is unacceptable) and the two sides may differ in size. Returns match with match[i] the receiver assigned to proposer i (or -1 if unmatched). The result is the proposer-optimal stable matching: it is stable (no blocking pair) and every proposer gets the best partner achievable in any stable matching.

Reference: Gale & Shapley, “College admissions and the stability of marriage”, Amer. Math. Monthly (1962).

Parameters:
Return type:

list[int]

tsp_held_karp(distance)[source]

Exact minimum-cost Hamiltonian cycle through all nodes (Held-Karp).

distance is an n x n matrix of arc costs (may be asymmetric); +inf marks a missing arc. NaN and -inf are refused (MXR-080-1903): the DP decides arc presence with np.isfinite, so both used to read as “missing”, silently answering for a smaller graph – and a -inf arc means an unboundedly profitable one, which is not the same fact as an absent one at all. Returns (cost, tour) where tour starts at node 0, visits every node once, and the cost includes the closing arc back to 0. Raises ValueError when the finite arcs admit no Hamiltonian cycle. The Held-Karp bitmask DP is exact in O(2^n n^2) time / O(2^n n) memory, so it is intended for small n (roughly <= 15-18); beyond that use a heuristic.

Parameters:

distance (Any)

Return type:

tuple[float, list[int]]

class Solution(value, objective)[source]

Bases: NamedTuple

One enumerated solution: value (the solution itself) and objective (its cost/score).

Parameters:
value: Any

Alias for field number 0

objective: float

Alias for field number 1

class SpanningTree(weights)[source]

Bases: Relation

Spanning trees of a weighted undirected graph, enumerated in increasing total edge weight.

The solution value is the list of (i, j) edges. Non-finite weights are forbidden edges.

Parameters:

weights (np.ndarray)

enumerator(k=None)[source]

Yield up to k spanning trees ordered by total edge weight.

Parameters:

k (int | None)

Return type:

Iterator[Solution]

class ViterbiPath(log_init, log_trans, log_obs)[source]

Bases: Relation

k most-likely hidden-state sequences of an HMM, enumerated in decreasing joint log-probability.

Standard Viterbi returns only the single best path; this reduces the trellis (nodes (t, s)) to a longest-log-prob path and yields the top k. The solution value is a length-T list of state indices. Delegates to mixle.enumeration.hmm_paths.hmm_best_paths(), whose backward Viterbi completion value is a tight admissible heuristic: emission log-densities may be positive (any density > 1), so a zero-heuristic sense="max" search would pop suboptimal goals first.

Parameters:

  • log_init (Any) – log p(state s at t=0), length S.

  • log_trans (Any) – log p(s' | s), shape (S, S).

  • log_obs (Any) – log p(observation_t | state s), shape (T, S) (emission log-likelihoods).

enumerator(k=None)[source]

Yield up to k hidden-state paths ordered by joint log probability.

Parameters:

k (int | None)

Return type:

Iterator[Solution]

best_first_paths(start, successors, is_goal=None, *, sense='min', heuristic=None, max_results=None, return_paths=True)[source]

Lazily enumerate goal states in monotone order of total additive cost/score.

Parameters:

  • start (Any) – The initial state.

  • successors (Callable[[Any], Iterable[tuple[Any, float]]]) – state -> iterable of (next_state, step) where step is the edge cost (sense="min") or edge score (sense="max").

  • is_goal (Callable[[Any], bool] | None) – Predicate; when true for a popped state it is emitted (and not expanded). None (the default) treats any sink – a state with no successors – as a goal, which covers DAG/trellis/edit-graph searches without a separate goal test.

  • sense (str) – "min" to minimise total cost (increasing order out) or "max" to maximise total score (decreasing order out).

  • heuristic (Callable[[Any], float] | None) – Optional admissible estimate of the remaining cost (a lower bound, min) or remaining score (an upper bound, max); 0 at goals. None is the always-admissible zero heuristic (uniform-cost search).

  • max_results (int | None) – Stop after yielding this many goals (None = exhaust the graph).

  • return_paths (bool) – Yield (path_list, total) when true, else (goal_state, total).

Yields:

(path_or_state, total) in best-first order; with an admissible heuristic the order is exact, and on a DAG/trellis the search enumerates k-best paths.

Return type:

Iterator[tuple[Any, float]]

nearest_first(start, neighbors, *, key=None, max_distance=None, max_results=None)[source]

Enumerate distinct states outward from start in increasing distance (Dijkstra).

The dual of best_first_paths(): instead of enumerating paths to goal states, this enumerates the reachable states themselves, each once, nearest first – an expanding metric “ball” around start. neighbors(state) -> iterable of (next_state, step_cost) with non-negative steps; key(state) gives a hashable identity for de-duplication (default: the state itself). The space may be infinite, so bound it with max_distance and/or max_results (or just consume the lazy iterator finitely).

Yields:

(state, distance) where distance is the shortest total cost from start, in nondecreasing order.

Parameters:
Return type:

Iterator[tuple[Any, float]]