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 and 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.

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]

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

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.

sense: str = 'min'
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]

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.

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(). 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 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).

sense: str = 'min'
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]

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

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).

sense: str = 'min'
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.

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 cheap. 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]

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

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.

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 cheapest arborescence in which every node is reachable from root; returns None if no such arborescence exists.

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.

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). Returns (cost, tour) where tour starts at node 0, visits every node once, and the cost includes the closing arc back to 0. 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)

sense: str = 'min'
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]

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.

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).

sense: str = 'max'
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]

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]]