From 59ef13e7d7d97beae3e11589a5b440f7c878f240 Mon Sep 17 00:00:00 2001 From: Berend <50019630+berendmarkhorst@users.noreply.github.com> Date: Fri, 4 Sep 2026 16:38:00 +0200 Subject: [PATCH] Harden exact solve status handling --- README.md | 1 + benchmarks/README.md | 14 + benchmarks/benchmark_dreyfus_memory.py | 196 ++++++++++++ docs/source/guide/solvers.md | 18 ++ docs/source/guide/variants.md | 8 + steinerpy/dreyfus_wagner.py | 86 +++-- steinerpy/dual_ascent.py | 10 + steinerpy/mathematical_model.py | 409 ++++++++++++++++++++---- steinerpy/objects.py | 147 ++++++++- steinerpy/pc_transform.py | 29 +- tests/test_dreyfus_wagner.py | 39 +++ tests/test_mwcsb.py | 1 + tests/test_nonnegative_costs.py | 127 ++++++++ tests/test_status_safety.py | 425 +++++++++++++++++++++++++ tests/test_steiner_problem.py | 7 +- 15 files changed, 1400 insertions(+), 117 deletions(-) create mode 100644 benchmarks/benchmark_dreyfus_memory.py create mode 100644 tests/test_nonnegative_costs.py create mode 100644 tests/test_status_safety.py diff --git a/README.md b/README.md index 9d5d8e3..d65871d 100644 --- a/README.md +++ b/README.md @@ -36,6 +36,7 @@ G.add_edge("B", "C", weight=2) G.add_edge("C", "D", weight=1) # One terminal group = Steiner tree; multiple groups = Steiner forest +# Edge/arc costs must be non-negative. solution = SteinerProblem(G, [["A", "D"]]).get_solution() print(f"Optimal cost: {solution.objective}") diff --git a/benchmarks/README.md b/benchmarks/README.md index e53f6e1..b76a0db 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -46,6 +46,20 @@ skipped and rows are appended, so a preempted sweep is restarted with the same command. Per-instance solve time is bounded by `--time-limit`; for hard wall-clock isolation on a cluster, the Slurm array runs one instance per task. +### Dreyfus-Wagner memory benchmark + +`benchmark_dreyfus_memory.py` compares the few-terminal dynamic program across +source checkouts with fixed generated graphs. Every repetition runs in a fresh +process so peak RSS is isolated, and the CSV records the Python, NetworkX, and +SciPy versions used: + +```bash +python benchmarks/benchmark_dreyfus_memory.py \ + --source-root /path/to/checkout --label candidate \ + --k 6,7,8,9,10 --sizes 300,1200 --densities 3,8 --repeats 3 \ + --output candidate.csv +``` + ## Output columns `instance, nodes, edges, terminals, opt, base_obj, base_rt, base_gap, diff --git a/benchmarks/benchmark_dreyfus_memory.py b/benchmarks/benchmark_dreyfus_memory.py new file mode 100644 index 0000000..7bf538f --- /dev/null +++ b/benchmarks/benchmark_dreyfus_memory.py @@ -0,0 +1,196 @@ +"""Reproducible Dreyfus-Wagner runtime and peak-RSS benchmark. + +Each repetition runs in a fresh subprocess so ``ru_maxrss`` is isolated. Use +the same interpreter with two source checkouts to compare implementations. +""" + +from __future__ import annotations + +import argparse +import csv +import json +import os +import resource +import statistics +import subprocess +import sys +import time +from pathlib import Path + + +DEFAULT_K = (6, 7, 8, 9, 10) +DEFAULT_SIZES = (300, 1200) +DEFAULT_DENSITIES = (3, 8) + + +def _connected_graph(nx, n: int, edge_factor: int, k: int, seed: int): + """Deterministic connected sparse graph with exactly about factor*n edges.""" + import random + + rng = random.Random(seed) + graph = nx.Graph() + graph.add_nodes_from(range(n)) + order = list(range(n)) + rng.shuffle(order) + for i in range(1, n): + parent = order[rng.randrange(i)] + graph.add_edge(order[i], parent, weight=rng.randint(1, 20)) + + target = min(n * (n - 1) // 2, edge_factor * n) + while graph.number_of_edges() < target: + u, v = rng.sample(range(n), 2) + if not graph.has_edge(u, v): + graph.add_edge(u, v, weight=rng.randint(1, 20)) + terminals = rng.sample(range(n), k) + return graph, terminals + + +def _worker(args) -> None: + sys.path.insert(0, str(Path(args.source_root).resolve())) + import networkx as nx + import scipy + + from steinerpy.dreyfus_wagner import dreyfus_wagner + + graph, terminals = _connected_graph( + nx, args.nodes, args.edge_factor, args.k, args.seed + ) + started = time.perf_counter() + objective, edges = dreyfus_wagner(graph, terminals) + runtime = time.perf_counter() - started + + selected = nx.Graph() + selected.add_nodes_from(terminals) + selected.add_edges_from(edges) + feasible = len(terminals) <= 1 or all( + nx.has_path(selected, terminals[0], terminal) for terminal in terminals[1:] + ) + raw_rss = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss + rss_bytes = raw_rss if sys.platform == "darwin" else raw_rss * 1024 + print( + json.dumps( + { + "runtime_s": runtime, + "peak_rss_mib": rss_bytes / (1024 * 1024), + "objective": objective, + "edge_count": len(edges), + "feasible": feasible, + "python": sys.version.split()[0], + "networkx": nx.__version__, + "scipy": scipy.__version__, + } + ) + ) + + +def _parse_ints(value: str): + return tuple(int(item) for item in value.split(",") if item) + + +def _coordinator(args) -> None: + rows = [] + script = str(Path(__file__).resolve()) + for k in args.k: + for n in args.sizes: + for factor in args.densities: + runtimes = [] + peaks = [] + reference = None + for _repetition in range(args.repeats): + seed = args.seed + 100_000 * k + 1_000 * n + factor + command = [ + sys.executable, + script, + "--worker", + "--source-root", + args.source_root, + "--k", + str(k), + "--nodes", + str(n), + "--edge-factor", + str(factor), + "--seed", + str(seed), + ] + completed = subprocess.run( + command, + check=True, + text=True, + capture_output=True, + env={**os.environ, "PYTHONHASHSEED": "0"}, + ) + result = json.loads(completed.stdout.strip().splitlines()[-1]) + if not result["feasible"]: + raise RuntimeError( + f"infeasible reconstruction for k={k}, n={n}, " + f"factor={factor}" + ) + signature = (result["objective"], result["edge_count"]) + if reference is None: + reference = signature + elif signature != reference: + raise RuntimeError( + f"non-deterministic result for k={k}, n={n}, " + f"factor={factor}: {signature} != {reference}" + ) + runtimes.append(result["runtime_s"]) + peaks.append(result["peak_rss_mib"]) + + assert reference is not None + rows.append( + { + "label": args.label, + "k": k, + "nodes": n, + "edge_factor": factor, + "edges": min(n * (n - 1) // 2, factor * n), + "seed": seed, + "repeats": args.repeats, + "median_runtime_s": statistics.median(runtimes), + "median_peak_rss_mib": statistics.median(peaks), + "objective": reference[0], + "python": result["python"], + "networkx": result["networkx"], + "scipy": result["scipy"], + } + ) + + fieldnames = list(rows[0]) + destination = Path(args.output) if args.output else None + handle = destination.open("w", newline="") if destination else sys.stdout + try: + writer = csv.DictWriter(handle, fieldnames=fieldnames) + writer.writeheader() + writer.writerows(rows) + finally: + if destination: + handle.close() + + +def _parser(): + parser = argparse.ArgumentParser() + parser.add_argument("--worker", action="store_true") + parser.add_argument("--source-root", default=".") + parser.add_argument("--label", default="current") + parser.add_argument("--output") + parser.add_argument("--repeats", type=int, default=3) + parser.add_argument("--k", type=_parse_ints, default=DEFAULT_K) + parser.add_argument("--sizes", type=_parse_ints, default=DEFAULT_SIZES) + parser.add_argument("--densities", type=_parse_ints, default=DEFAULT_DENSITIES) + parser.add_argument("--nodes", type=int) + parser.add_argument("--edge-factor", type=int) + parser.add_argument("--seed", type=int, default=20260904) + return parser + + +if __name__ == "__main__": + args = _parser().parse_args() + if args.worker: + if isinstance(args.k, tuple): + if len(args.k) != 1: + raise SystemExit("--worker requires one --k") + args.k = args.k[0] + _worker(args) + else: + _coordinator(args) diff --git a/docs/source/guide/solvers.md b/docs/source/guide/solvers.md index a6f82e8..fc9e3c1 100644 --- a/docs/source/guide/solvers.md +++ b/docs/source/guide/solvers.md @@ -19,6 +19,24 @@ solution = SteinerProblem(graph, terminal_groups).get_solution(solver="gurobi") Both solvers implement the same cut-based (DO-D) formulation from Markhorst et al. (2025) and produce identical optimal solutions. Gurobi may be faster on larger instances because callbacks avoid repeated re-solves from scratch. +## Input and certificate guarantees + +All edge-cost Steiner variants require non-negative values in the selected +edge/arc weight attribute. This is required by the shortest-path reductions +and exact formulations. +Negative node weights remain valid for maximum-weight connected-subgraph +problems. Inputs with a negative edge or arc cost raise `ValueError` at +construction time instead of risking a disconnected negative-cost cycle in the +reported solution. + +A time limit is not an optimality certificate. If a solve stops early with an +independently validated feasible incumbent, `get_solution()` may return it +with a nonzero or unknown gap (`math.inf`). If no valid incumbent exists—or +the last incumbent predates newly added connectivity cuts—the values are +discarded and the public solve raises `RuntimeError`. In particular, +`solution.gap == 0` is reported only after the solver or an independent exact +bound proves optimality. + ## Enumerating multiple optimal solutions `get_solution()` returns exactly one optimal Steiner tree. When multiple diff --git a/docs/source/guide/variants.md b/docs/source/guide/variants.md index 246aa83..41f2e92 100644 --- a/docs/source/guide/variants.md +++ b/docs/source/guide/variants.md @@ -20,6 +20,14 @@ Simply pass a list of terminal lists as `terminal_groups` — one list for a tre Several of these variants implement the "further related problems" of Chapter 5 of D. Rehfeldt's PhD thesis (*Faster algorithms for Steiner tree and related problems*, TU Berlin 2021), reusing the same directed-cut kernel by transformation. +All edge and arc costs used by these exact Steiner formulations must be +non-negative. This restriction does not apply to node *weights* in +`MaxWeightConnectedSubgraph` or +`BudgetedMaxWeightConnectedSubgraph`, where negative values mean that a node +is an undesirable connector and are part of the problem semantics. MWCSPB also +treats edge attributes as topology-only; its objective and budget use node +weights and node costs. + ```python from steinerpy import ( PartialTerminalSteinerProblem, FullTerminalSteinerProblem, GroupSteinerProblem, diff --git a/steinerpy/dreyfus_wagner.py b/steinerpy/dreyfus_wagner.py index 3e1d5d7..693225e 100644 --- a/steinerpy/dreyfus_wagner.py +++ b/steinerpy/dreyfus_wagner.py @@ -14,6 +14,11 @@ *virtual source* connected to every vertex with its merged label as arc cost — a standard trick to run Dijkstra with initial potentials. scipy is required (callers gate on :data:`steinerpy._fastgraph.HAS_SCIPY`). + +Only the subset label arrays are retained during the forward pass. +Reconstruction recomputes split choices and shortest-path predecessors for the +states actually visited by the final tree, trading a small amount of final +work for lower peak memory. """ from __future__ import annotations @@ -60,6 +65,17 @@ def dreyfus_wagner(graph, terminals: List[Hashable], weight: str = "weight" if not HAS_SCIPY: raise RuntimeError("dreyfus_wagner requires scipy") + negative = [ + (u, v, data.get(weight, 1)) + for u, v, data in graph.edges(data=True) + if data.get(weight, 1) < 0 + ] + if negative: + raise ValueError( + "dreyfus_wagner requires non-negative edge weights; " + f"found {negative[:5]}{' ...' if len(negative) > 5 else ''}." + ) + terminals = list(dict.fromkeys(terminals)) if len(terminals) <= 1: return 0.0, [] @@ -88,25 +104,19 @@ def dreyfus_wagner(graph, terminals: List[Hashable], weight: str = "weight" kb = len(base) full = (1 << kb) - 1 - # Per subset S (index = bitmask over `base`): - # labels[S][v] = l(S, v), length n+1 (the virtual entry is unused); - # grow_pred[S] = predecessor array of the grow Dijkstra — `virtual` - # marks the base vertex where the merged subtree sits - # (for singletons: predecessors of the plain Dijkstra); - # split[S][v] = the canonical submask chosen by the merge at v. + # Forward DP stores only l(S, v). Predecessors and split choices are needed + # solely along the final tree, so reconstruction recomputes those few states + # instead of retaining two additional subset-indexed array families. labels = [None] * (full + 1) - grow_pred = [None] * (full + 1) - split = [None] * (full + 1) graph_csr = csr_matrix((base_costs, (base_tails, base_heads)), shape=(n + 1, n + 1)) # Base case: l({t}, v) = d(t, v) — one plain Dijkstra per base terminal. for i, ti in enumerate(base): - dist, pred = _sp_dijkstra(graph_csr, directed=True, indices=ti, - return_predecessors=True) - labels[1 << i] = dist - grow_pred[1 << i] = pred + labels[1 << i] = _sp_dijkstra( + graph_csr, directed=True, indices=ti + ) # Every grow step uses the same sparsity pattern: the fixed graph arcs plus # one virtual-source arc to each real vertex. Build that pattern once and @@ -132,50 +142,74 @@ def dreyfus_wagner(graph, terminals: List[Hashable], weight: str = "weight" # Merge: m(S, v) = min over canonical splits S' (containing the lowest # set bit, so each unordered pair is tried once) of l(S') + l(S \ S'). best = np.full(n + 1, np.inf) - choice = np.zeros(n + 1, dtype=np.int64) sub = (mask - 1) & mask while sub: if sub & low: cand = labels[sub] + labels[mask ^ sub] upd = cand < best best[upd] = cand[upd] - choice[upd] = sub sub = (sub - 1) & mask # Grow: l(S, v) = min_u m(S, u) + d(u, v) — a Dijkstra from a virtual # source whose arc to u costs m(S, u). grow_csr.data[virtual_slice] = best[virtual_heads] - dist, pred = _sp_dijkstra(grow_csr, directed=True, indices=virtual, - return_predecessors=True) - labels[mask] = dist - grow_pred[mask] = pred - split[mask] = choice + labels[mask] = _sp_dijkstra( + grow_csr, directed=True, indices=virtual + ) total = labels[full][root] if not np.isfinite(total): return float("inf"), [] - # Reconstruction: walk grow predecessors back to the subtree base, then - # expand the merge recorded there into its two sub-subsets. + # Reconstruction: recompute a predecessor array only for states on the + # chosen tree. Walk grow predecessors back to the subtree base, then + # recompute the canonical split at that one vertex. edge_keys = set() stack = [(full, root)] while stack: mask, v = stack.pop() if mask & (mask - 1) == 0: source = base[mask.bit_length() - 1] + _dist, pred = _sp_dijkstra( + graph_csr, directed=True, indices=source, + return_predecessors=True, + ) w = v while w != source: - u = int(grow_pred[mask][w]) + u = int(pred[w]) edge_keys.add((u, w) if u < w else (w, u)) w = u continue + + low = mask & -mask + best = np.full(n + 1, np.inf) + sub = (mask - 1) & mask + while sub: + if sub & low: + best = np.minimum(best, labels[sub] + labels[mask ^ sub]) + sub = (sub - 1) & mask + + grow_csr.data[virtual_slice] = best[virtual_heads] + _dist, pred = _sp_dijkstra( + grow_csr, directed=True, indices=virtual, + return_predecessors=True, + ) w = v while True: - u = int(grow_pred[mask][w]) + u = int(pred[w]) if u == virtual: - sub = int(split[mask][w]) - stack.append((sub, w)) - stack.append((mask ^ sub, w)) + chosen = 0 + chosen_cost = np.inf + sub = (mask - 1) & mask + while sub: + if sub & low: + candidate = labels[sub][w] + labels[mask ^ sub][w] + if candidate < chosen_cost: + chosen_cost = candidate + chosen = sub + sub = (sub - 1) & mask + stack.append((chosen, w)) + stack.append((mask ^ chosen, w)) break edge_keys.add((u, w) if u < w else (w, u)) w = u diff --git a/steinerpy/dual_ascent.py b/steinerpy/dual_ascent.py index 752a529..5e7db92 100644 --- a/steinerpy/dual_ascent.py +++ b/steinerpy/dual_ascent.py @@ -534,6 +534,16 @@ def dual_ascent(steiner_problem, weight: Optional[str] = None) -> DualAscentResu """Run dual ascent + primal heuristic on ``steiner_problem.graph``.""" weight = weight or steiner_problem.weight graph = steiner_problem.graph + negative = [ + (u, v, data.get(weight, 1)) + for u, v, data in graph.edges(data=True) + if data.get(weight, 1) < 0 + ] + if negative: + raise ValueError( + "dual_ascent requires non-negative edge/arc costs; " + f"found {negative[:5]}{' ...' if len(negative) > 5 else ''}." + ) arcs = list(steiner_problem.arcs) edges = list(steiner_problem.edges) groups = steiner_problem.terminal_groups diff --git a/steinerpy/mathematical_model.py b/steinerpy/mathematical_model.py index 552b324..51f6f4d 100644 --- a/steinerpy/mathematical_model.py +++ b/steinerpy/mathematical_model.py @@ -1,6 +1,5 @@ import highspy as hp import logging -import math import networkx as nx import os import time @@ -33,6 +32,50 @@ def _resolve_threads(threads) -> int: return 0 +def _highs_outcome(model: hp.HighsModel) -> Tuple[str, bool]: + """Classify a HiGHS termination without reading an invalid solution. + + The boolean reports whether a primal incumbent is available. In particular, + a time-limit status is *incomplete* even when HiGHS happened to close the + numerical MIP gap; only ``kOptimal`` is treated as an optimality proof. + """ + status = model.getModelStatus() + if status == hp.HighsModelStatus.kInfeasible: + return "infeasible", False + if status in ( + hp.HighsModelStatus.kUnbounded, + hp.HighsModelStatus.kUnboundedOrInfeasible, + hp.HighsModelStatus.kNotset, + hp.HighsModelStatus.kLoadError, + hp.HighsModelStatus.kModelError, + hp.HighsModelStatus.kPresolveError, + hp.HighsModelStatus.kSolveError, + hp.HighsModelStatus.kPostsolveError, + hp.HighsModelStatus.kUnknown, + ): + return "incomplete", False + try: + has_incumbent = bool(model.getSolution().value_valid) + except Exception: # pragma: no cover - defensive for old highspy releases + has_incumbent = False + if not has_incumbent: + return "incomplete", False + if status == hp.HighsModelStatus.kOptimal: + return "optimal", True + return "incomplete", True + + +def _gurobi_outcome(model, GRB) -> Tuple[str, bool]: + """Gurobi counterpart of :func:`_highs_outcome`.""" + if model.SolCount <= 0: + if model.Status == GRB.INFEASIBLE: + return "infeasible", False + return "incomplete", False + if model.Status == GRB.OPTIMAL: + return "optimal", True + return "incomplete", True + + # Separation parallelism: number of worker threads for the per-terminal min-cut # computations. scipy's maximum_flow releases the GIL, so threads overlap. def _sep_thread_count() -> int: @@ -155,6 +198,32 @@ def _incident_edges(edges) -> Dict: return incident +def _selection_graph(steiner_problem, selected_edges): + graph_type = ( + nx.DiGraph + if isinstance(steiner_problem.graph, nx.DiGraph) + else nx.Graph + ) + selected = graph_type() + selected.add_nodes_from(steiner_problem.nodes) + selected.add_edges_from(selected_edges) + return selected + + +def _groups_reachable(steiner_problem, selected_edges, skipped=None) -> bool: + """Validate required terminal reachability in an extracted incumbent.""" + skipped = skipped or set() + selected = _selection_graph(steiner_problem, selected_edges) + for group_id, group in enumerate(steiner_problem.terminal_groups): + root = steiner_problem.roots[group_id] + for terminal in dict.fromkeys(group): + if terminal == root or (group_id, terminal) in skipped: + continue + if not nx.has_path(selected, root, terminal): + return False + return True + + def get_terminals(terminal_group: List[List]) -> List: """ Turns a nested list of terminals into a list of terminals. @@ -911,6 +980,7 @@ def _ret(gap, runtime, objective, selected_edges, status): # Phase 2 — the integer cut loop. converged = False + feasible_incumbent = False # ever_solved distinguishes "the time budget ran out before a single MIP # solve completed" (e.g. time_limit=0) from a real incumbent: without it, # falling through to read model.variableValue()/getObjectiveValue() below @@ -930,15 +1000,12 @@ def _ret(gap, runtime, objective, selected_edges, status): # separate cuts from — reading the empty solution would yield spurious # "violations" and loop forever. Stop the cut loop instead of hanging. status = model.getModelStatus() - if status == hp.HighsModelStatus.kInfeasible: + outcome, has_incumbent = _highs_outcome(model) + if outcome == "infeasible": runtime = time.time() - start_time logging.info("Cut loop: model proved infeasible.") return _ret(float("inf"), runtime, float("inf"), [], "infeasible") - if status in ( - hp.HighsModelStatus.kObjectiveBound, - hp.HighsModelStatus.kUnbounded, - hp.HighsModelStatus.kUnboundedOrInfeasible, - ) or not model.getSolution().value_valid: + if not has_incumbent: runtime = time.time() - start_time logging.warning( "Cut loop stopped early: model status %s with no valid primal " @@ -954,6 +1021,7 @@ def _ret(gap, runtime, objective, selected_edges, status): # with no violated cuts; treating it as converged would falsely report # gap 0 (the false-optimal bug fixed in solve_sap_highs). converged = status == hp.HighsModelStatus.kOptimal + feasible_incumbent = True break # feasible w.r.t. all cut constraints # Add each violated cut as a new constraint: sum(y2[k,a] for a in cut) >= z[k,l] @@ -966,8 +1034,11 @@ def _ret(gap, runtime, objective, selected_edges, status): runtime = time.time() - start_time logging.info(f"Runtime: {runtime:.2f} seconds") - if not ever_solved: - logging.warning("Cut loop: time limit exhausted before any solve completed.") + if not ever_solved or not feasible_incumbent: + logging.warning( + "Cut loop: no connectivity-valid incumbent was available when the " + "solve stopped." + ) return _ret(float("inf"), runtime, float("inf"), [], "incomplete") selected_edges = [e for e in steiner_problem.edges if model.variableValue(x[e]) > 0.5] @@ -976,9 +1047,9 @@ def _ret(gap, runtime, objective, selected_edges, status): gap = model.getInfo().mip_gap result_status = "optimal" else: - # Global time limit hit before all connectivity cuts were separated: - # selected_edges may be disconnected and the relaxation objective is only a - # lower bound, so do not report a (spurious) ~0 MIP gap. + # The incumbent passed independent cut separation but optimality was not + # proved. Do not surface a solver-reported zero gap for a non-optimal + # termination. gap = float("inf") result_status = "incomplete" @@ -1088,8 +1159,14 @@ def add_prize_collecting_constraints(model: hp.HighsModel, steiner_problem: 'Pri model.addConstr(total_penalties <= steiner_problem.penalty_budget) -def run_prize_collecting_model(model: hp.HighsModel, steiner_problem: 'PrizeCollectingProblem', - x: hp.HighsVarType, node_vars: Dict, penalty_vars: Dict) -> Tuple[float, float, float, List[Tuple], List[str], Dict]: +def run_prize_collecting_model( + model: hp.HighsModel, + steiner_problem: 'PrizeCollectingProblem', + x: hp.HighsVarType, + node_vars: Dict, + penalty_vars: Dict, + return_status: bool = False, +) -> Tuple: """ Solve prize collecting model and extract solution. """ @@ -1110,26 +1187,58 @@ def run_prize_collecting_model(model: hp.HighsModel, steiner_problem: 'PrizeColl for penalty_var in penalty_vars.values(): objective_expr += penalty_var * penalty_cost + def _ret(gap, runtime, objective, selected_edges, selected_nodes, penalties, + status): + result = (gap, runtime, objective, selected_edges, selected_nodes, penalties) + return result + (status,) if return_status else result + # Minimize the objective model.minimize(objective_expr) - - logging.info(f"Runtime: {model.getRunTime():.2f} seconds") + + runtime = model.getRunTime() + logging.info(f"Runtime: {runtime:.2f} seconds") + + status, has_incumbent = _highs_outcome(model) + if not has_incumbent: + return _ret(float("inf"), runtime, float("inf"), [], [], {}, status) # Extract solution selected_edges = [e for e in steiner_problem.edges if model.variableValue(x[e]) > 0.5] selected_nodes = [node for node in steiner_problem.nodes if model.variableValue(node_vars[node]) > 0.5] penalties = {} + penalized_terminals = set() for (group_id, terminal), var in penalty_vars.items(): var_value = model.variableValue(var) if var_value > 0.5: + penalized_terminals.add((group_id, terminal)) penalties[f"group_{group_id}_{terminal}"] = penalty_cost * var_value + + selected_graph = _selection_graph(steiner_problem, selected_edges) + roots = set(steiner_problem.roots) + nodes_connected = all( + node in roots or any( + nx.has_path(selected_graph, root, node) for root in roots + ) + for node in selected_nodes + ) + candidate_valid = _groups_reachable( + steiner_problem, selected_edges, penalized_terminals + ) + candidate_valid = candidate_valid and nodes_connected + if not candidate_valid: + logging.warning( + "Prize-collecting solve returned an invalid incumbent; discarding it." + ) + return _ret( + float("inf"), runtime, float("inf"), [], [], {}, "incomplete" + ) - gap = model.getInfo().mip_gap - runtime = model.getRunTime() + gap = model.getInfo().mip_gap if status == "optimal" else float("inf") objective = model.getObjectiveValue() - - return gap, runtime, objective, selected_edges, selected_nodes, penalties + + return _ret(gap, runtime, objective, selected_edges, selected_nodes, + penalties, status) def build_budget_model(steiner_problem: 'BaseSteinerProblem', time_limit: float = 300, logfile: str = "", threads=None) -> Tuple: @@ -1194,7 +1303,8 @@ def build_budget_model(steiner_problem: 'BaseSteinerProblem', time_limit: float def run_budget_model(model: hp.HighsModel, steiner_problem: 'BaseSteinerProblem', - x: Dict, penalty_vars: Dict) -> Tuple[float, float, int, List[Tuple], Dict]: + x: Dict, penalty_vars: Dict, + return_status: bool = False) -> Tuple: """ Solve budget-constrained model: minimize number of unconnected terminals. @@ -1206,24 +1316,64 @@ def run_budget_model(model: hp.HighsModel, steiner_problem: 'BaseSteinerProblem' """ logging.info("Started running the budget-constrained model...") + def _ret(gap, runtime, connected_count, selected_edges, penalties, status): + result = (gap, runtime, connected_count, selected_edges, penalties) + return result + (status,) if return_status else result + + # A budget instance with only group roots has a constant-zero objective and + # needs no solve (passing a Python 0 to highspy.minimize crashes). + if not penalty_vars: + total_terminals = sum(len(g) for g in steiner_problem.terminal_groups) + return _ret(0.0, 0.0, total_terminals, [], {}, "optimal") + model.minimize(sum(penalty_vars.values())) - logging.info(f"Runtime: {model.getRunTime():.2f} seconds") + runtime = model.getRunTime() + logging.info(f"Runtime: {runtime:.2f} seconds") + + status, has_incumbent = _highs_outcome(model) + if not has_incumbent: + return _ret(float("inf"), runtime, 0, [], {}, status) selected_edges = [e for e in steiner_problem.edges if model.variableValue(x[e]) > 0.5] penalties = {} + penalized_terminals = set() for (group_id, terminal), var in penalty_vars.items(): if model.variableValue(var) > 0.5: + penalized_terminals.add((group_id, terminal)) penalties[f"group_{group_id}_{terminal}"] = 1 + edge_cost = sum( + steiner_problem.graph.edges[e][steiner_problem.weight] + for e in selected_edges + ) + degree_ok = True + if getattr(steiner_problem, "max_degree", None) is not None: + degree_ok = all( + degree <= steiner_problem.max_degree + for _node, degree in _selection_graph( + steiner_problem, selected_edges + ).degree() + ) + candidate_valid = edge_cost <= steiner_problem.budget + 1e-7 + candidate_valid = candidate_valid and degree_ok + candidate_valid = candidate_valid and _groups_reachable( + steiner_problem, selected_edges, penalized_terminals + ) + if not candidate_valid: + logging.warning( + "Budget-constrained solve returned an invalid incumbent; " + "discarding it." + ) + return _ret(float("inf"), runtime, 0, [], {}, "incomplete") + total_terminals = sum(len(g) for g in steiner_problem.terminal_groups) connected_count = total_terminals - len(penalties) - gap = model.getInfo().mip_gap - runtime = model.getRunTime() + gap = model.getInfo().mip_gap if status == "optimal" else float("inf") - return gap, runtime, connected_count, selected_edges, penalties + return _ret(gap, runtime, connected_count, selected_edges, penalties, status) # --------------------------------------------------------------------------- @@ -1512,19 +1662,35 @@ def _cut_callback(cb_model, where): logging.info(f"Gurobi runtime: {runtime:.2f} seconds") - if model.SolCount == 0: + outcome, has_incumbent = _gurobi_outcome(model, GRB) + if not has_incumbent: # No feasible solution found. Distinguish a proven-infeasible model # from one that simply ran out of time before finding any incumbent # (e.g. time_limit=0) -- both previously looked identical (an # infinite objective), so a probe that merely timed out with no # incumbent was silently treated as proof there is no solution. - if model.Status in (GRB.INFEASIBLE, GRB.INF_OR_UNBD): + if outcome == "infeasible": return _ret(float("inf"), runtime, float("inf"), [], "infeasible") return _ret(float("inf"), runtime, float("inf"), [], "incomplete") selected_edges = [e for e in steiner_problem.edges if x[e].X > 0.5] objective = model.ObjVal - if model.Status == GRB.OPTIMAL: + # A callback can be interrupted at the deadline. Independently separate + # the final incumbent before exposing it; SolCount alone does not certify + # that every lazy connectivity constraint was processed. + y2_vals = { + (group_id, a): y2[(group_id, a)].X + for group_id in group_indices for a in steiner_problem.arcs + } + z_vals = {key: var.X for key, var in z.items()} + if find_violated_cuts_from_values(steiner_problem, y2_vals, z_vals): + logging.warning( + "Gurobi stopped with an incumbent that violates connectivity; " + "discarding it." + ) + return _ret(float("inf"), runtime, float("inf"), [], "incomplete") + + if outcome == "optimal": gap = model.MIPGap status = "optimal" else: @@ -1564,8 +1730,7 @@ def _sap_indegree_into(view) -> Dict: def solve_sap_highs(view, time_limit: float = 300, logfile: str = "", fixing=None, da_cuts=None, da_ub=None, primal=None, - threads=None - ) -> Tuple[float, float, float, List[Tuple]]: + threads=None, return_status: bool = False) -> Tuple: """Solve a single-group SAP by HiGHS + iterative directed-cut generation. ``view`` is any object exposing ``graph`` (DiGraph), ``arcs``, ``nodes``, @@ -1580,6 +1745,10 @@ def solve_sap_highs(view, time_limit: float = 300, logfile: str = "", :returns: ``(gap, runtime, objective, selected_arcs)``. """ + def _ret(gap, runtime, objective, selected, status): + result = (gap, runtime, objective, selected) + return result + (status,) if return_status else result + model = make_model(time_limit, logfile, threads=threads) arcs = list(view.arcs) root = view.roots[0] @@ -1612,7 +1781,7 @@ def solve_sap_highs(view, time_limit: float = 300, logfile: str = "", # as a termination target inside the cut loop: a loose dual-ascent UB makes a # re-solve stop kOptimal at a feasible-but-suboptimal incumbent, which the loop # then reports as "proven optimal" (false-optimal bug observed on PCSPG P400). - # da_ub is still used below only to report an honest gap on a non-proven solve. + # da_ub remains in the signature for backend symmetry and compatibility. obj = sum(xa[a] * view.graph.edges[a][view.weight] for a in arcs) @@ -1661,12 +1830,8 @@ def solve_sap_highs(view, time_limit: float = 300, logfile: str = "", # Phase 2 — the integer cut loop. converged = False - _STOP = ( - hp.HighsModelStatus.kInfeasible, - hp.HighsModelStatus.kObjectiveBound, - hp.HighsModelStatus.kUnbounded, - hp.HighsModelStatus.kUnboundedOrInfeasible, - ) + feasible_incumbent = False + outcome = "incomplete" while True: # Bound the *total* cut-generation time, not just each individual re-solve. # The loop can add many rounds of Steiner cuts, and each model.minimize() @@ -1680,7 +1845,8 @@ def solve_sap_highs(view, time_limit: float = 300, logfile: str = "", model.setOptionValue("time_limit", remaining) model.minimize(obj) status = model.getModelStatus() - if status in _STOP or not model.getSolution().value_valid: + outcome, has_incumbent = _highs_outcome(model) + if not has_incumbent: break col_value = model.getSolution().col_value xa_vals = {(0, a): col_value[xa[a].index] for a in arcs} @@ -1691,37 +1857,72 @@ def solve_sap_highs(view, time_limit: float = 300, logfile: str = "", # suboptimal incumbent with no violated cuts; treating that as converged # would stamp a non-optimal tree "proven optimal" (observed on P400). converged = status == hp.HighsModelStatus.kOptimal + feasible_incumbent = True break + # The current solution is optimal only for the still-incomplete + # relaxation. Once a violated cut is found it is not a feasible SAP + # incumbent, even if the deadline expires before the next re-solve. + outcome = "incomplete" for (_k, _l, cut_arcs) in violated: if cut_arcs: model.addConstr(sum(xa[a] for a in cut_arcs) >= 1) runtime = time.time() - start - selected = [a for a in arcs if model.variableValue(xa[a]) > 0.5] - objective = model.getObjectiveValue() + if feasible_incumbent: + selected = [a for a in arcs if model.variableValue(xa[a]) > 0.5] + objective = model.getObjectiveValue() + else: + # The last solver values may pre-date connectivity cuts that were just + # added. Never expose them. A supplied dual-ascent primal is an + # independently feasible fallback, so validate and use it when possible. + selected = list(primal or []) + if selected and not _sap_candidate_feasible(view, selected): + selected = [] + if selected: + objective = sum(view.graph.edges[a][view.weight] for a in selected) + outcome = "incomplete" + else: + return _ret(float("inf"), runtime, float("inf"), [], outcome) + if converged: gap = model.getInfo().mip_gap + result_status = "optimal" else: - # Time limit hit before all connectivity cuts were separated: the model is - # a relaxation whose optimum is only a lower bound, so its MIP gap would be - # a spurious ~0. Report an honest gap against the dual-ascent feasible upper - # bound instead (the caller maps the best valid component of `selected`). - if da_ub is not None and math.isfinite(da_ub) and da_ub > 0: - gap = max(0.0, (da_ub - objective) / max(1.0, abs(da_ub))) - else: - gap = float("inf") - return gap, runtime, objective, selected + # The incumbent is connectivity-valid, but the solver did not prove + # optimality. Its MIP gap (including a numerical zero) is not an exact + # certificate for the still-incomplete cut model. + gap = float("inf") + result_status = "incomplete" + return _ret(gap, runtime, objective, selected, result_status) + + +def _sap_candidate_feasible(view, selected) -> bool: + """Validate SAP connectivity and the arborescence indegree constraints.""" + chosen = set(selected) + root = view.roots[0] + indegree: Dict = {} + for u, v in chosen: + if (u, v) not in view.arcs or v == root: + return False + indegree[v] = indegree.get(v, 0) + 1 + if indegree[v] > 1: + return False + vals = {(0, a): 1.0 if a in chosen else 0.0 for a in view.arcs} + return not find_violated_cuts_from_values(view, vals, {(0, 0): 1.0}) def solve_sap_gurobi(view, time_limit: float = 300, logfile: str = "", fixing=None, da_cuts=None, da_ub=None, primal=None, - threads=None - ) -> Tuple[float, float, float, List[Tuple]]: + threads=None, return_status: bool = False) -> Tuple: """Gurobi branch-and-cut counterpart of :func:`solve_sap_highs`. Connectivity is separated lazily inside a callback, mirroring :func:`run_model_gurobi`. """ + def _ret(gap, runtime, objective, selected, status): + result = (gap, runtime, objective, selected) + return result + (status,) if return_status else result + _check_gurobipy() import gurobipy as gp from gurobipy import GRB @@ -1794,10 +1995,21 @@ def _cb(cb_model, where): model.optimize(_cb) runtime = time.time() - start - if model.SolCount == 0: - return float("inf"), runtime, float("inf"), [] - selected = [a for a in arcs if xa[a].X > 0.5] - return model.MIPGap, runtime, model.ObjVal, selected + outcome, has_incumbent = _gurobi_outcome(model, GRB) + selected = [a for a in arcs if xa[a].X > 0.5] if has_incumbent else [] + candidate_valid = has_incumbent and _sap_candidate_feasible(view, selected) + if has_incumbent and not candidate_valid: + selected = [] + outcome = "incomplete" + if not candidate_valid and primal and _sap_candidate_feasible(view, primal): + selected = list(primal) + outcome = "incomplete" + candidate_valid = True + if not candidate_valid: + return _ret(float("inf"), runtime, float("inf"), [], outcome) + objective = sum(view.graph.edges[a][view.weight] for a in selected) + gap = model.MIPGap if outcome == "optimal" else float("inf") + return _ret(gap, runtime, objective, selected, outcome) # --------------------------------------------------------------------------- @@ -1824,6 +2036,24 @@ def _undirected_from_arcs(used_arcs: List[Tuple]) -> List[Tuple]: return edges +def _mwcsb_candidate_feasible( + steiner_problem, selected_edges, selected_nodes +) -> bool: + chosen = set(selected_nodes) + root = steiner_problem.roots[0] + if root not in chosen: + return False + if any(u not in chosen or v not in chosen for u, v in selected_edges): + return False + selected = nx.Graph() + selected.add_nodes_from(chosen) + selected.add_edges_from(selected_edges) + if any(not nx.has_path(selected, root, node) for node in chosen): + return False + spent = sum(steiner_problem.node_costs.get(v, 0) for v in chosen) + return spent <= steiner_problem.node_budget + 1e-7 + + def build_mwcsb_model(steiner_problem, time_limit: float = 300, logfile: str = "", threads=None) -> Tuple: """Build the HiGHS MWCSPB model (see module note above). @@ -1866,7 +2096,8 @@ def build_mwcsb_model(steiner_problem, time_limit: float = 300, logfile: str = " return model, x, y1, y2, z, node_vars -def run_mwcsb_model(model, steiner_problem, y1: Dict, node_vars: Dict) -> Tuple: +def run_mwcsb_model(model, steiner_problem, y1: Dict, node_vars: Dict, + return_status: bool = False) -> Tuple: """Solve the HiGHS MWCSPB model and extract the connected subgraph. :return: ``(gap, runtime, mwcs_weight, selected_edges, selected_nodes)`` where @@ -1877,20 +2108,40 @@ def run_mwcsb_model(model, steiner_problem, y1: Dict, node_vars: Dict) -> Tuple: objective_expr = sum(-nw.get(v, 0.0) * node_vars[v] for v in steiner_problem.nodes) model.minimize(objective_expr) - status = model.getModelStatus() - if status in ( - hp.HighsModelStatus.kInfeasible, - hp.HighsModelStatus.kUnbounded, - hp.HighsModelStatus.kUnboundedOrInfeasible, - ) or not model.getSolution().value_valid: - return float("inf"), model.getRunTime(), float("-inf"), [], [] + def _ret(gap, runtime, objective, selected_edges, selected_nodes, status): + result = (gap, runtime, objective, selected_edges, selected_nodes) + return result + (status,) if return_status else result + + status, has_incumbent = _highs_outcome(model) + if not has_incumbent: + return _ret( + float("inf"), model.getRunTime(), float("-inf"), [], [], status + ) selected_nodes = [v for v in steiner_problem.nodes if model.variableValue(node_vars[v]) > 0.5] - used_arcs = [a for a in steiner_problem.arcs if model.variableValue(y1[a]) > 0.5] + chosen = set(selected_nodes) + # The root has no indegree constraint. An otherwise unused arc entering it may + # therefore be one in a degenerate optimum, even though its tail is not selected. + # Such an arc is not part of the node-induced MWCSPB solution. + used_arcs = [ + a + for a in steiner_problem.arcs + if a[0] in chosen and a[1] in chosen + and model.variableValue(y1[a]) > 0.5 # noqa: W503 + ] selected_edges = _undirected_from_arcs(used_arcs) mwcs_weight = sum(nw.get(v, 0.0) for v in selected_nodes) + if not _mwcsb_candidate_feasible( + steiner_problem, selected_edges, selected_nodes + ): + return _ret( + float("inf"), model.getRunTime(), float("-inf"), [], [], + "incomplete", + ) - return model.getInfo().mip_gap, model.getRunTime(), mwcs_weight, selected_edges, selected_nodes + gap = model.getInfo().mip_gap if status == "optimal" else float("inf") + return _ret(gap, model.getRunTime(), mwcs_weight, selected_edges, + selected_nodes, status) def build_mwcsb_model_gurobi(steiner_problem, time_limit: float = 300, logfile: str = "", @@ -1956,7 +2207,8 @@ def build_mwcsb_model_gurobi(steiner_problem, time_limit: float = 300, logfile: return model, x, y1, y2, z, node_vars -def run_mwcsb_model_gurobi(model, steiner_problem, y1: Dict, node_vars: Dict) -> Tuple: +def run_mwcsb_model_gurobi(model, steiner_problem, y1: Dict, node_vars: Dict, + return_status: bool = False) -> Tuple: """Gurobi counterpart of :func:`run_mwcsb_model`.""" _check_gurobipy() import gurobipy as gp @@ -1971,12 +2223,29 @@ def run_mwcsb_model_gurobi(model, steiner_problem, y1: Dict, node_vars: Dict) -> model.optimize() runtime = time.time() - start - if model.SolCount == 0: - return float("inf"), runtime, float("-inf"), [], [] + def _ret(gap, runtime, objective, selected_edges, selected_nodes, status): + result = (gap, runtime, objective, selected_edges, selected_nodes) + return result + (status,) if return_status else result + + status, has_incumbent = _gurobi_outcome(model, GRB) + if not has_incumbent: + return _ret(float("inf"), runtime, float("-inf"), [], [], status) selected_nodes = [v for v in steiner_problem.nodes if node_vars[v].X > 0.5] - used_arcs = [a for a in steiner_problem.arcs if y1[a].X > 0.5] + chosen = set(selected_nodes) + used_arcs = [ + a + for a in steiner_problem.arcs + if a[0] in chosen and a[1] in chosen and y1[a].X > 0.5 + ] selected_edges = _undirected_from_arcs(used_arcs) mwcs_weight = sum(nw.get(v, 0.0) for v in selected_nodes) + if not _mwcsb_candidate_feasible( + steiner_problem, selected_edges, selected_nodes + ): + return _ret( + float("inf"), runtime, float("-inf"), [], [], "incomplete" + ) - return model.MIPGap, runtime, mwcs_weight, selected_edges, selected_nodes + gap = model.MIPGap if status == "optimal" else float("inf") + return _ret(gap, runtime, mwcs_weight, selected_edges, selected_nodes, status) diff --git a/steinerpy/objects.py b/steinerpy/objects.py index f5e917a..9f765b7 100644 --- a/steinerpy/objects.py +++ b/steinerpy/objects.py @@ -14,6 +14,23 @@ logger = logging.getLogger(__name__) +def _validate_nonnegative_edge_costs(graph: nx.Graph, weight: str, problem: str) -> None: + """Reject costs that invalidate shortest-path and exact-model assumptions.""" + negative = [ + (u, v, data.get(weight, 1)) + for u, v, data in graph.edges(data=True) + if data.get(weight, 1) < 0 + ] + if negative: + preview = negative[:5] + suffix = " ..." if len(negative) > len(preview) else "" + raise ValueError( + f"{problem} requires non-negative edge/arc costs in attribute " + f"{weight!r}; found {preview}{suffix}. Negative node weights remain " + "supported by maximum-weight connected-subgraph variants." + ) + + def node_split_graph( graph: nx.Graph, terminal_groups: List[List], @@ -53,6 +70,8 @@ def node_split_graph( class BaseSteinerProblem: + _requires_nonnegative_edge_costs = True + def __init__(self, graph: nx.Graph, terminal_groups: List[List], weight="weight", preprocess=True, **kwargs): """ Initialize the SteinerProblem (can be tree or forest). @@ -61,6 +80,14 @@ def __init__(self, graph: nx.Graph, terminal_groups: List[List], weight="weight" :param terminal_groups: nested list of terminals. :param weight: edge attribute specified by this string as the edge weight. """ + # ``enumeration_safe`` has an older, stricter positive-cost contract and + # its own public error message. Let that check below retain precedence. + if self._requires_nonnegative_edge_costs and not ( + preprocess and kwargs.get("enumeration_safe", False) + ): + _validate_nonnegative_edge_costs( + graph, weight, type(self).__name__ + ) self.original_graph = graph self.preprocess = preprocess # Opt-in dual-ascent bound reduction (removes provably non-optimal edges @@ -203,6 +230,22 @@ def _da_eligible(self) -> bool: return False return True + def _connects_terminal_groups(self, edges: List[Tuple]) -> bool: + """Independently validate the public connectivity contract.""" + graph_type = nx.DiGraph if isinstance(self.graph, nx.DiGraph) else nx.Graph + selected = graph_type() + selected.add_nodes_from(self.nodes) + selected.add_edges_from(edges) + for group in self.terminal_groups: + terminals = list(dict.fromkeys(group)) + if len(terminals) <= 1: + continue + root = terminals[0] + if any(not nx.has_path(selected, root, terminal) + for terminal in terminals[1:]): + return False + return True + def _dw_eligible(self) -> bool: """Whether the Dreyfus–Wagner dynamic program can solve this instance. @@ -230,7 +273,8 @@ def _dw_eligible(self) -> bool: k = len(set(self.terminal_groups[0])) if k < 2 or k > dw_max_terminals(): return False - # ~3 arrays of 2^(k-1) x (n+1) doubles/ints; keep the footprint modest. + # One label array per subset, plus bounded reconstruction temporaries; + # keep the exponential footprint modest. if (1 << (k - 1)) * (self.graph.number_of_nodes() + 1) > 50_000_000: return False return True @@ -606,9 +650,15 @@ def get_solution(self, time_limit: float = 300, log_file: str = "", solver: str model, x, y1, y2, z, f, penalty_vars = build_budget_model( self, time_limit=time_limit, logfile=log_file, threads=threads ) - gap, runtime, connected_count, selected_edges, penalties = run_budget_model( - model, self, x, penalty_vars + gap, runtime, connected_count, selected_edges, penalties, _status = run_budget_model( + model, self, x, penalty_vars, return_status=True ) + if _status != "optimal" and connected_count == 0: + reason = ( + "was proved infeasible" if _status == "infeasible" + else "stopped before finding a valid feasible incumbent" + ) + raise RuntimeError(f"Budget-constrained solve {reason}.") if self.preprocess: original_selected_edges = map_solution_to_original( @@ -677,14 +727,15 @@ def get_solution(self, time_limit: float = 300, log_file: str = "", solver: str # Optional dual-ascent accelerator: lower bound + primal heuristic + # reduced-cost variable fixing. Early-exits when proven optimal. + import math as _math use_da = self.dual_ascent if dual_ascent is None else dual_ascent fixing = None + da = None da_primal = None da_cuts = None da_ub = None if use_da and self._da_eligible(): import time as _time - import math as _math from .dual_ascent import ( dual_ascent as _run_da, reduced_cost_fixing, steiner_cuts, ) @@ -711,7 +762,7 @@ def get_solution(self, time_limit: float = 300, log_file: str = "", solver: str apply_fixes_gurobi(model, x, y1, y2, fixing) set_gurobi_cutoff(model, da_ub) gap, runtime, objective, selected_edges, _status = run_model_gurobi(model, self, x, y2, z, return_status=True) - if da_cuts is not None and _math.isinf(objective): + if da_cuts is not None and _status == "infeasible": # The dual-ascent acceleration (cutoff / fixing / seeded cuts) # over-constrained the model into infeasibility although the # instance is feasible; re-solve from a clean model. @@ -734,12 +785,39 @@ def get_solution(self, time_limit: float = 300, log_file: str = "", solver: str _m, _x, _p = model, x, da_primal reapply_start = lambda: set_highs_warm_start(_m, _x, _p) # noqa: E731 gap, runtime, objective, selected_edges, _status = run_model(model, self, x, y2, z, reapply_start=reapply_start, return_status=True) - if da_cuts is not None and _math.isinf(objective): + if da_cuts is not None and _status == "infeasible": # The acceleration over-constrained a feasible instance into # infeasibility; re-solve from a clean, un-accelerated model. model, x, y1, y2, z = build_model(self, time_limit=time_limit, logfile=log_file, threads=threads) gap, runtime, objective, selected_edges, _status = run_model(model, self, x, y2, z, return_status=True) + if _status == "incomplete" and not _math.isfinite(objective): + if da is None or not da.feasible or not _math.isfinite(da.upper_bound): + raise RuntimeError( + "Exact solve stopped before finding a connectivity-valid " + "feasible incumbent." + ) + # Dual ascent supplies an independently feasible incumbent and lower + # bound. Use it rather than returning stale values from a cut round. + selected_edges = list(da.primal_edges) + objective = da.upper_bound + if not _math.isfinite(da.lower_bound): + gap = float("inf") + else: + gap = max( + 0.0, + (da.upper_bound - da.lower_bound) + / max(1.0, abs(da.upper_bound)), + ) + + if _math.isfinite(objective) and not self._connects_terminal_groups( + selected_edges + ): + raise RuntimeError( + "Exact solver returned an incumbent that does not connect every " + "terminal group; the result was discarded." + ) + # Map solution back to original graph if preprocessing was used if self.preprocess: original_selected_edges = map_solution_to_original(selected_edges, self.reduction_tracker, self.graph) @@ -1106,6 +1184,7 @@ def __init__(self, graph: nx.Graph, terminal_groups: List[List], partial_termina weight: str = "weight", **kwargs): if isinstance(graph, nx.DiGraph): raise ValueError("PartialTerminalSteinerProblem requires an undirected graph.") + _validate_nonnegative_edge_costs(graph, weight, type(self).__name__) self.partial_terminals = set(partial_terminals) all_terminals = {t for group in terminal_groups for t in group} missing = self.partial_terminals - all_terminals @@ -1576,10 +1655,10 @@ def _pc_exact_solution(self, time_limit, log_file, solver, threads=None) -> 'Pri da_ub = da.upper_bound solve = solve_sap_gurobi if solver == "gurobi" else solve_sap_highs - gap, _rt, _sap_obj, sap_arcs = solve( + gap, _rt, _sap_obj, sap_arcs, _status = solve( view, time_limit=time_limit, logfile=log_file, fixing=fixing, da_cuts=da_cuts, da_ub=da_ub, primal=da_primal, - threads=threads, + threads=threads, return_status=True, ) edges, nodes, pcstp_obj = map_sap_solution_to_pcstp(ctx, sap_arcs) @@ -1689,9 +1768,16 @@ def get_solution(self, time_limit: float = 300, log_file: str = "", self, time_limit=time_limit, logfile=log_file, threads=threads ) - gap, runtime, objective, selected_edges, selected_nodes, penalties = run_prize_collecting_model( - model, self, x, node_vars, penalty_vars + (gap, runtime, objective, selected_edges, selected_nodes, penalties, + _status) = run_prize_collecting_model( + model, self, x, node_vars, penalty_vars, return_status=True ) + if _status != "optimal" and objective == float("inf"): + reason = ( + "was proved infeasible" if _status == "infeasible" + else "stopped before finding a valid feasible incumbent" + ) + raise RuntimeError(f"Prize-collecting solve {reason}.") # Map solution back to original graph if preprocessing was used if self.preprocess: @@ -1955,6 +2041,14 @@ def __init__(self, graph: nx.Graph, terminal_groups: List[List], node_weights: D Note: graph preprocessing is not supported for node-weighted problems because the node-splitting transformation produces a directed graph internally. """ + negative_nodes = [(v, cost) for v, cost in node_weights.items() if cost < 0] + if negative_nodes: + raise ValueError( + "NodeWeightedSteinerProblem requires non-negative node costs; " + f"found {negative_nodes[:5]}" + f"{' ...' if len(negative_nodes) > 5 else ''}. Use " + "MaxWeightConnectedSubgraph for semantically negative node weights." + ) self.node_weights = node_weights self.original_terminal_groups_nw = terminal_groups @@ -1989,6 +2083,17 @@ def get_solution(self, time_limit: float = 300, log_file: str = "", solver: str model, x, y1, y2, z = build_model(self, time_limit=time_limit, logfile=log_file, threads=threads) gap, runtime, objective, _, _status = run_model(model, self, x, y2, z, return_status=True) + if _status == "incomplete" and objective == float("inf"): + raise RuntimeError( + "Node-weighted exact solve stopped before finding a " + "connectivity-valid feasible incumbent." + ) + if objective == float("inf"): + return NodeWeightedSolution( + gap=gap, runtime=runtime, objective=objective, + selected_edges=[], original_selected_edges=[], selected_nodes=[], + ) + # Use arc (y1) variables for the actual directed tree structure instead of edge # (x) variables, to avoid degenerate zero-weight cross-edges being included. # Also exclude arcs pointing INTO root nodes (roots have no parents in an arborescence). @@ -2184,6 +2289,10 @@ class BudgetedMaxWeightConnectedSubgraph(MaxWeightConnectedSubgraph): problems: From theory to practice*, PhD thesis, TU Berlin, Ch. 5.6. """ + # Original edge attributes are topology-only in MWCSPB: neither the + # node-weight objective nor the node-cost budget reads them. + _requires_nonnegative_edge_costs = False + def __init__(self, graph: nx.Graph, node_weights: Dict, node_costs: Dict, node_budget: float, root=None, weight: str = "weight", **kwargs): """ @@ -2213,16 +2322,25 @@ def get_solution(self, time_limit: float = 300, log_file: str = "", model, x, y1, y2, z, node_vars = build_mwcsb_model_gurobi( self, time_limit=time_limit, logfile=log_file, threads=threads ) - gap, runtime, mwcs_weight, selected_edges, selected_nodes = run_mwcsb_model_gurobi( - model, self, y1, node_vars + (gap, runtime, mwcs_weight, selected_edges, selected_nodes, + _status) = run_mwcsb_model_gurobi( + model, self, y1, node_vars, return_status=True ) else: model, x, y1, y2, z, node_vars = build_mwcsb_model( self, time_limit=time_limit, logfile=log_file, threads=threads ) - gap, runtime, mwcs_weight, selected_edges, selected_nodes = run_mwcsb_model( - model, self, y1, node_vars + (gap, runtime, mwcs_weight, selected_edges, selected_nodes, + _status) = run_mwcsb_model( + model, self, y1, node_vars, return_status=True + ) + + if _status != "optimal" and mwcs_weight == float("-inf"): + reason = ( + "was proved infeasible" if _status == "infeasible" + else "stopped before finding a valid feasible incumbent" ) + raise RuntimeError(f"Budgeted MWCS solve {reason}.") total_prize = sum(max(0.0, self._mwcs_node_weights.get(v, 0.0)) for v in selected_nodes) return PrizeCollectingSolution( @@ -2405,6 +2523,7 @@ def __init__(self, graph: nx.DiGraph, root, terminals: List, hop_limit: int, """ if not isinstance(graph, nx.DiGraph): raise ValueError("HopConstrainedSteinerProblem requires a directed graph (nx.DiGraph).") + _validate_nonnegative_edge_costs(graph, weight, type(self).__name__) # Drop outgoing arcs of every non-root terminal (thesis: delta^+_S(t) = 0). non_root_terminals = {t for t in terminals if t != root} diff --git a/steinerpy/pc_transform.py b/steinerpy/pc_transform.py index c137429..5831ac5 100644 --- a/steinerpy/pc_transform.py +++ b/steinerpy/pc_transform.py @@ -55,6 +55,19 @@ def _term_label(t) -> Arc: return ("__pc_term__", t) +def _validate_nonnegative_costs(graph, weight: str, transform: str) -> None: + negative = [ + (u, v, data.get(weight, 1)) + for u, v, data in graph.edges(data=True) + if data.get(weight, 1) < 0 + ] + if negative: + raise ValueError( + f"{transform} requires non-negative edge/arc costs; found " + f"{negative[:5]}{' ...' if len(negative) > 5 else ''}." + ) + + # --------------------------------------------------------------------------- # Transform context # --------------------------------------------------------------------------- @@ -137,6 +150,7 @@ def transform_pcstp_to_sap( :returns: a :class:`PCTransform`. """ + _validate_nonnegative_costs(graph, weight, "transform_pcstp_to_sap") proper: Set = set() nonproper: Set = set() for v in graph.nodes(): @@ -241,6 +255,9 @@ def transform_directed_pcstp_to_sap( raise ValueError( "transform_directed_pcstp_to_sap requires a directed graph (nx.DiGraph)." ) + _validate_nonnegative_costs( + graph, weight, "transform_directed_pcstp_to_sap" + ) if root is not None and root not in graph: raise ValueError(f"root {root!r} is not in the graph.") @@ -314,8 +331,9 @@ def transform_mwcsp_to_pcstp( ) -> Tuple[nx.Graph, Dict, float]: """MWCSP ``(graph, node_weights)`` -> equivalent classic PCSTP. - Rehfeldt & Koch (2020), Sec. 2.2: let ``w0 = min_v w(v)``. Define edge costs - ``c(e) := -w0`` (``> 0`` when some weight is negative) and prizes + Rehfeldt & Koch (2020), Sec. 2.2: let + ``w0 = min(0, min_v w(v))``. Define edge costs ``c(e) := -w0`` + (``> 0`` when some weight is negative) and prizes ``p(v) := w(v) - w0 >= 0``. Maximising ``sum_{v in S} w(v)`` over connected ``S`` is then equivalent to minimising the PCSTP objective, and the original weight is recovered by ``mwcsp_weight = mwcsp_const - pcstp_obj`` with @@ -326,12 +344,15 @@ def transform_mwcsp_to_pcstp( """ if not node_weights: raise ValueError("node_weights must be non-empty for MWCSP.") - w0 = min(node_weights.get(v, 0) for v in graph.nodes()) + # Clamping at zero matters when all node weights are positive. The algebraic + # tree transformation still works with a positive shift, but it would create + # negative edge costs and invalidate the SAP shortest-path/cut machinery. + w0 = min(0, min(node_weights.get(v, 0) for v in graph.nodes())) n = graph.number_of_nodes() pc_graph = nx.Graph() pc_graph.add_nodes_from(graph.nodes()) - edge_cost = -w0 # >= 0; 0 only in the degenerate all-nonnegative case + edge_cost = -w0 # >= 0; zero when all node weights are non-negative for u, v in graph.edges(): pc_graph.add_edge(u, v, **{weight: edge_cost}) diff --git a/tests/test_dreyfus_wagner.py b/tests/test_dreyfus_wagner.py index f39fb79..7e97bb1 100644 --- a/tests/test_dreyfus_wagner.py +++ b/tests/test_dreyfus_wagner.py @@ -80,6 +80,20 @@ def test_zero_weight_edges(): assert _spans(edges, [0, 2]) +def test_tied_shortest_paths_reconstruct_a_valid_optimum(): + # Two equal routes between each side of the terminal set exercise both + # Dijkstra predecessor ties and tied merge choices during reconstruction. + G = nx.cycle_graph(6) + for u, v in G.edges(): + G[u][v]["weight"] = 1 + + cost, edges = dreyfus_wagner(G, [0, 2, 4]) + + assert cost == pytest.approx(4.0) + assert _edge_cost(G, edges) == pytest.approx(cost) + assert _spans(edges, [0, 2, 4]) + + def test_disconnected_terminals(): G = nx.Graph() G.add_edge(0, 1, weight=1) @@ -107,6 +121,31 @@ def counting_csr_matrix(*args, **kwargs): assert len(calls) == 2 # base graph plus reusable virtual-source graph +def test_forward_dp_does_not_retain_predecessor_families(monkeypatch): + """Predecessors are requested only for states visited in reconstruction.""" + module = importlib.import_module("steinerpy.dreyfus_wagner") + real_dijkstra = module._sp_dijkstra + predecessor_calls = 0 + + def counting_dijkstra(*args, **kwargs): + nonlocal predecessor_calls + if kwargs.get("return_predecessors", False): + predecessor_calls += 1 + return real_dijkstra(*args, **kwargs) + + monkeypatch.setattr(module, "_sp_dijkstra", counting_dijkstra) + graph, terminals = _random_instance(40, 120, 6, seed=405) + + cost, edges = module.dreyfus_wagner(graph, terminals) + + assert math.isfinite(cost) + assert _spans(edges, terminals) + # The old implementation requested predecessors for all 2^(k-1)-1 states + # (31 here). Label-only reconstruction visits at most a binary tree of + # 2*(k-1)-1 states. + assert predecessor_calls <= 2 * (len(terminals) - 1) - 1 + + def test_disconnected_grow_step_with_infinite_virtual_weights(): """The reusable virtual row may contain inf for unreachable vertices.""" graph = nx.Graph() diff --git a/tests/test_mwcsb.py b/tests/test_mwcsb.py index 19a60f1..ebec163 100644 --- a/tests/test_mwcsb.py +++ b/tests/test_mwcsb.py @@ -73,6 +73,7 @@ def test_zero_budget_keeps_only_root(solver): ).get_solution(solver=solver) assert sol.objective == pytest.approx(5.0) assert set(sol.selected_nodes) == {"r"} + assert sol.edges == [] def test_highs_and_gurobi_match(): diff --git a/tests/test_nonnegative_costs.py b/tests/test_nonnegative_costs.py new file mode 100644 index 0000000..91cdcf8 --- /dev/null +++ b/tests/test_nonnegative_costs.py @@ -0,0 +1,127 @@ +"""Non-negative edge-cost preconditions for exact Steiner algorithms.""" + +from types import SimpleNamespace + +import networkx as nx +import pytest + +from steinerpy import ( + BudgetedMaxWeightConnectedSubgraph, + DirectedSteinerProblem, + HopConstrainedSteinerProblem, + MaxWeightConnectedSubgraph, + NodeWeightedSteinerProblem, + PartialTerminalSteinerProblem, + PrizeCollectingProblem, + SteinerProblem, +) +from steinerpy.dreyfus_wagner import dreyfus_wagner +from steinerpy.dual_ascent import dual_ascent +from steinerpy.pc_transform import ( + transform_directed_pcstp_to_sap, + transform_mwcsp_to_pcstp, + transform_pcstp_to_sap, +) + + +@pytest.mark.parametrize("preprocess", [False, True]) +def test_disconnected_negative_cycle_is_rejected(preprocess): + """A disconnected negative cycle must never lower a Steiner objective.""" + graph = nx.Graph() + graph.add_edge("s", "m", weight=4) + graph.add_edge("m", "t", weight=4) + graph.add_edge("a", "b", weight=-3) + graph.add_edge("b", "c", weight=-3) + graph.add_edge("c", "a", weight=-3) + + with pytest.raises(ValueError, match="non-negative edge/arc"): + SteinerProblem(graph, [["s", "t"]], preprocess=preprocess) + + +def test_directed_and_prize_collecting_negative_edges_are_rejected(): + directed = nx.DiGraph() + directed.add_edge("r", "t", cost=-1) + with pytest.raises(ValueError, match="attribute 'cost'"): + DirectedSteinerProblem(directed, root="r", terminals=["t"], weight="cost") + + graph = nx.Graph() + graph.add_edge(0, 1, weight=-1) + with pytest.raises(ValueError, match="non-negative edge/arc"): + PrizeCollectingProblem(graph, [[0]], {0: 1, 1: 1}, penalty_cost=0) + + +def test_transformed_variants_validate_before_dropping_or_changing_edges(): + graph = nx.Graph() + graph.add_edge(0, 1, weight=-1) + graph.add_edge(1, 2, weight=2) + with pytest.raises(ValueError, match="non-negative"): + PartialTerminalSteinerProblem(graph, [[0, 1, 2]], partial_terminals=[0]) + + directed = nx.DiGraph() + directed.add_edge("r", "t", weight=-1) + with pytest.raises(ValueError, match="non-negative"): + HopConstrainedSteinerProblem(directed, root="r", terminals=["t"], hop_limit=1) + + +def test_direct_algorithm_and_transform_entry_points_reject_negative_costs(): + graph = nx.Graph() + graph.add_edge(0, 1, weight=-1) + with pytest.raises(ValueError, match="non-negative"): + dreyfus_wagner(graph, [0, 1]) + with pytest.raises(ValueError, match="non-negative"): + transform_pcstp_to_sap(graph, {0: 1, 1: 1}) + + view = SimpleNamespace(graph=graph, weight="weight") + with pytest.raises(ValueError, match="non-negative"): + dual_ascent(view) + + directed = nx.DiGraph() + directed.add_edge("r", "t", weight=-1) + with pytest.raises(ValueError, match="non-negative"): + transform_directed_pcstp_to_sap(directed, {"r": 0, "t": 1}, root="r") + + +def test_zero_and_missing_edge_costs_remain_supported(): + graph = nx.Graph() + graph.add_edge(0, 1, weight=0) + graph.add_edge(1, 2) # missing costs keep the historical default of 1 + + problem = SteinerProblem(graph, [[0, 2]], preprocess=False) + + assert problem.graph[0][1]["weight"] == 0 + assert "weight" not in problem.graph[1][2] + + +def test_negative_node_weights_remain_supported_for_mwcs_variants(): + graph = nx.path_graph(3) + nx.set_edge_attributes(graph, 0, "weight") + weights = {0: 5, 1: -2, 2: 4} + + mwcs = MaxWeightConnectedSubgraph(graph, weights, root=0) + topology_only = graph.copy() + nx.set_edge_attributes(topology_only, -99, "weight") + budgeted = BudgetedMaxWeightConnectedSubgraph( + topology_only, + weights, + {0: 0, 1: 1, 2: 1}, + node_budget=2, + root=0, + ) + + assert mwcs._mwcs_node_weights[1] == -2 + assert budgeted._mwcs_node_weights[1] == -2 + + +def test_all_positive_mwcs_transform_keeps_edge_costs_nonnegative(): + graph = nx.path_graph(3) + transformed, prizes, constant = transform_mwcsp_to_pcstp(graph, {0: 5, 1: 2, 2: 4}) + + assert set(nx.get_edge_attributes(transformed, "weight").values()) == {0} + assert prizes == {0: 5, 1: 2, 2: 4} + assert constant == 11 + + +def test_negative_node_costs_are_rejected_for_node_weighted_steiner(): + graph = nx.path_graph(3) + with pytest.raises(ValueError, match="non-negative node costs"): + NodeWeightedSteinerProblem(graph, [[0, 2]], node_weights={0: 0, 1: -1, 2: 0}) diff --git a/tests/test_status_safety.py b/tests/test_status_safety.py new file mode 100644 index 0000000..e5c3fdf --- /dev/null +++ b/tests/test_status_safety.py @@ -0,0 +1,425 @@ +"""Regression tests for time-limited and otherwise incomplete exact solves.""" + +import importlib.util +import math +from types import SimpleNamespace + +import highspy as hp +import networkx as nx +import pytest + +import steinerpy.mathematical_model as mm +import steinerpy.objects as objects +from steinerpy import ( + BudgetedMaxWeightConnectedSubgraph, + NodeWeightedSteinerProblem, + PrizeCollectingProblem, + SteinerProblem, +) + + +def _cycle_problem(n=20): + graph = nx.cycle_graph(n) + nx.set_edge_attributes(graph, 1, "weight") + return SteinerProblem(graph, [[0, n // 4, n // 2]], preprocess=False) + + +def _require_gurobi(): + if importlib.util.find_spec("gurobipy") is None: + pytest.skip("gurobipy is not installed") + try: + import gurobipy as gp + + env = gp.Env(empty=True) + env.setParam("OutputFlag", 0) + env.start() + model = gp.Model(env=env) + model.dispose() + env.dispose() + except Exception: + pytest.skip("Gurobi license is not available") + + +def test_tiny_positive_limit_has_no_bogus_core_incumbent(monkeypatch): + """A positive limit can expire before the first MIP solve starts.""" + monkeypatch.setenv("STEINERPY_LP_CUT_ROUNDS", "0") + problem = _cycle_problem() + model, x, _y1, y2, z = mm.build_model(problem, time_limit=1e-9) + + gap, _runtime, objective, edges, status = mm.run_model( + model, problem, x, y2, z, return_status=True + ) + + assert status == "incomplete" + assert math.isinf(gap) and math.isinf(objective) + assert edges == [] + + +def test_deadline_after_violated_cut_discards_previous_incumbent(monkeypatch): + """The pre-cut incumbent is invalid once its violated cut has been added.""" + monkeypatch.setenv("STEINERPY_LP_CUT_ROUNDS", "0") + problem = _cycle_problem(8) + model, x, _y1, y2, z = mm.build_model(problem, time_limit=1.0) + + real_separate = mm.find_violated_cuts + separated = [] + + def recording_separation(*args, **kwargs): + cuts = real_separate(*args, **kwargs) + separated.append(cuts) + return cuts + + # start=0; permit the first solve at t=0; expire before the re-solve at t=2. + times = iter((0.0, 0.0, 2.0, 2.0)) + monkeypatch.setattr(mm, "time", SimpleNamespace(time=lambda: next(times))) + monkeypatch.setattr(mm, "find_violated_cuts", recording_separation) + + gap, _runtime, objective, edges, status = mm.run_model( + model, problem, x, y2, z, return_status=True + ) + + assert separated and separated[0] + assert status == "incomplete" + assert math.isinf(gap) and math.isinf(objective) + assert edges == [] + + +def test_public_get_solution_rejects_disconnected_incomplete_result( + monkeypatch, +): + graph = nx.path_graph(3) + nx.set_edge_attributes(graph, 1, "weight") + problem = SteinerProblem(graph, [[0, 2]], preprocess=False) + monkeypatch.setenv("STEINERPY_DW_MAX_TERMINALS", "0") + monkeypatch.setattr( + objects, + "run_model", + lambda *args, **kwargs: (math.inf, 0.01, 0.0, [], "incomplete"), + ) + + with pytest.raises(RuntimeError, match="does not connect"): + problem.get_solution(decompose=False) + + +def test_public_get_solution_accepts_valid_unproven_incumbent(monkeypatch): + graph = nx.path_graph(3) + nx.set_edge_attributes(graph, 1, "weight") + problem = SteinerProblem(graph, [[0, 2]], preprocess=False) + monkeypatch.setenv("STEINERPY_DW_MAX_TERMINALS", "0") + monkeypatch.setattr( + objects, + "run_model", + lambda *args, **kwargs: ( + math.inf, + 0.01, + 2.0, + [(0, 1), (1, 2)], + "incomplete", + ), + ) + + solution = problem.get_solution(decompose=False) + + assert solution.edges == [(0, 1), (1, 2)] + assert solution.objective == pytest.approx(2.0) + assert math.isinf(solution.gap) + + +def test_preprocessed_tiny_limit_reports_no_valid_solution(monkeypatch): + """Preprocessing/back-mapping must not turn a missing incumbent into a tree.""" + graph = nx.complete_graph(14) + nx.set_edge_attributes(graph, 1, "weight") + monkeypatch.setenv("STEINERPY_DW_MAX_TERMINALS", "0") + monkeypatch.setenv("STEINERPY_LP_CUT_ROUNDS", "0") + problem = SteinerProblem( + graph, + [[0, 4, 9]], + preprocess=True, + heavy=False, + contract_terminals=False, + bound_based=False, + ) + + with pytest.raises(RuntimeError, match="before finding"): + problem.get_solution(time_limit=1e-9, decompose=False) + + +def test_core_infeasible_status_is_distinct(monkeypatch): + monkeypatch.setenv("STEINERPY_LP_CUT_ROUNDS", "0") + graph = nx.Graph() + graph.add_edge(0, 1, weight=1) + graph.add_edge(2, 3, weight=1) + problem = SteinerProblem(graph, [[0, 3]], preprocess=False) + model, x, _y1, y2, z = mm.build_model(problem, time_limit=5.0) + + gap, _runtime, objective, edges, status = mm.run_model( + model, problem, x, y2, z, return_status=True + ) + + assert status == "infeasible" + assert math.isinf(gap) and math.isinf(objective) + assert edges == [] + + +def test_specialized_highs_runners_reject_no_incumbent(): + graph = nx.cycle_graph(20) + nx.set_edge_attributes(graph, 1, "weight") + + pc = PrizeCollectingProblem(graph, [[0, 10]], {v: 1 for v in graph}, penalty_cost=1) + model, x, _y1, _y2, _z, _f, node_vars, penalty_vars = ( + mm.build_prize_collecting_model(pc, time_limit=1e-9) + ) + pc_result = mm.run_prize_collecting_model( + model, pc, x, node_vars, penalty_vars, return_status=True + ) + assert pc_result[-1] == "incomplete" + assert math.isinf(pc_result[0]) and math.isinf(pc_result[2]) + assert pc_result[3:6] == ([], [], {}) + + budget = SteinerProblem(graph, [[0, 10]], preprocess=False, budget=5) + model, x, _y1, _y2, _z, _f, penalty_vars = mm.build_budget_model( + budget, time_limit=1e-9 + ) + budget_result = mm.run_budget_model( + model, budget, x, penalty_vars, return_status=True + ) + assert budget_result[-1] == "incomplete" + assert math.isinf(budget_result[0]) + assert budget_result[2:5] == (0, [], {}) + + mwcsb = BudgetedMaxWeightConnectedSubgraph( + graph, + {v: 1 for v in graph}, + {v: 1 for v in graph}, + node_budget=10, + root=0, + ) + model, _x, y1, _y2, _z, node_vars = mm.build_mwcsb_model(mwcsb, time_limit=1e-9) + mwcsb_result = mm.run_mwcsb_model(model, mwcsb, y1, node_vars, return_status=True) + assert mwcsb_result[-1] == "incomplete" + assert math.isinf(mwcsb_result[0]) and mwcsb_result[2] == -math.inf + assert mwcsb_result[3:5] == ([], []) + + +def test_specialized_infeasible_and_valid_unproven_statuses(monkeypatch): + graph = nx.path_graph(3) + nx.set_edge_attributes(graph, 1, "weight") + + # Root cost alone exceeds the budget: the full-flow MWCSPB model is + # genuinely infeasible, rather than merely lacking an incumbent. + infeasible = BudgetedMaxWeightConnectedSubgraph( + graph, + {0: 5, 1: 1, 2: 1}, + {0: 2, 1: 1, 2: 1}, + node_budget=1, + root=0, + ) + model, _x, y1, _y2, _z, node_vars = mm.build_mwcsb_model(infeasible, time_limit=5) + result = mm.run_mwcsb_model(model, infeasible, y1, node_vars, return_status=True) + assert result[-1] == "infeasible" + assert result[2] == -math.inf + + # Solve a feasible specialized model, then force only the reported solver + # termination to kTimeLimit. Values remain a valid incumbent, but gap=0 is + # forbidden because optimality is no longer reported as proved. + feasible = BudgetedMaxWeightConnectedSubgraph( + graph, + {0: 5, 1: 1, 2: 10}, + {0: 0, 1: 1, 2: 1}, + node_budget=2, + root=0, + ) + model, _x, y1, _y2, _z, node_vars = mm.build_mwcsb_model(feasible, time_limit=5) + model.minimize( + sum( + -feasible._mwcs_node_weights.get(v, 0) * node_vars[v] + for v in feasible.nodes + ) + ) + assert model.getSolution().value_valid + monkeypatch.setattr(model, "getModelStatus", lambda: hp.HighsModelStatus.kTimeLimit) + result = mm.run_mwcsb_model(model, feasible, y1, node_vars, return_status=True) + assert result[-1] == "incomplete" + assert math.isinf(result[0]) + assert set(result[4]) == {0, 1, 2} + + +def test_prize_and_budget_valid_unproven_incumbents(monkeypatch): + graph = nx.path_graph(3) + nx.set_edge_attributes(graph, 1, "weight") + + pc = PrizeCollectingProblem(graph, [[0, 2]], {0: 2, 1: 1, 2: 2}, penalty_cost=10) + model, x, _y1, _y2, _z, _f, node_vars, penalty_vars = ( + mm.build_prize_collecting_model(pc, time_limit=5) + ) + first = mm.run_prize_collecting_model( + model, pc, x, node_vars, penalty_vars, return_status=True + ) + assert first[-1] == "optimal" + monkeypatch.setattr(model, "getModelStatus", lambda: hp.HighsModelStatus.kTimeLimit) + result = mm.run_prize_collecting_model( + model, pc, x, node_vars, penalty_vars, return_status=True + ) + assert result[-1] == "incomplete" + assert math.isinf(result[0]) and math.isfinite(result[2]) + assert result[3] + + budget = SteinerProblem(graph, [[0, 2]], preprocess=False, budget=2) + model, x, _y1, _y2, _z, _f, penalty_vars = mm.build_budget_model( + budget, time_limit=5 + ) + first = mm.run_budget_model(model, budget, x, penalty_vars, return_status=True) + assert first[-1] == "optimal" + monkeypatch.setattr(model, "getModelStatus", lambda: hp.HighsModelStatus.kTimeLimit) + result = mm.run_budget_model(model, budget, x, penalty_vars, return_status=True) + assert result[-1] == "incomplete" + assert math.isinf(result[0]) + assert result[2] == 2 and result[3] + + +def test_sap_tiny_positive_limit_and_valid_primal_fallback(monkeypatch): + monkeypatch.setenv("STEINERPY_LP_CUT_ROUNDS", "0") + graph = nx.DiGraph() + graph.add_edge("r", "a", weight=1) + graph.add_edge("a", "t", weight=1) + view = objects.DirectedSteinerProblem(graph, root="r", terminals=["t"]) + + no_incumbent = mm.solve_sap_highs(view, time_limit=1e-9, return_status=True) + assert no_incumbent[-1] == "incomplete" + assert math.isinf(no_incumbent[0]) and math.isinf(no_incumbent[2]) + assert no_incumbent[3] == [] + + fallback = mm.solve_sap_highs( + view, + time_limit=1e-9, + primal=[("r", "a"), ("a", "t")], + return_status=True, + ) + assert fallback[-1] == "incomplete" + assert math.isinf(fallback[0]) + assert fallback[2] == pytest.approx(2.0) + assert fallback[3] == [("r", "a"), ("a", "t")] + + +def test_sap_deadline_after_violated_cut_discards_incumbent(monkeypatch): + monkeypatch.setenv("STEINERPY_LP_CUT_ROUNDS", "0") + graph = nx.DiGraph() + graph.add_edge("r", "a", weight=1) + graph.add_edge("a", "t", weight=1) + view = objects.DirectedSteinerProblem(graph, root="r", terminals=["t"]) + times = iter((0.0, 0.0, 2.0, 2.0)) + monkeypatch.setattr(mm, "time", SimpleNamespace(time=lambda: next(times))) + + result = mm.solve_sap_highs(view, time_limit=1.0, return_status=True) + + assert result[-1] == "incomplete" + assert math.isinf(result[0]) and math.isinf(result[2]) + assert result[3] == [] + + +def test_public_specialized_paths_reject_missing_incumbents(monkeypatch): + graph = nx.cycle_graph(30) + nx.set_edge_attributes(graph, 1, "weight") + + pc = PrizeCollectingProblem(graph, [[0, 15]], {v: 1 for v in graph}, penalty_cost=1) + with pytest.raises(RuntimeError, match="before finding"): + pc.get_solution(time_limit=1e-9) + + budget = SteinerProblem(graph, [[0, 15]], preprocess=False, budget=5) + with pytest.raises(RuntimeError, match="before finding"): + budget.get_solution(time_limit=1e-9) + + mwcsb = BudgetedMaxWeightConnectedSubgraph( + graph, + {v: 1 for v in graph}, + {v: 1 for v in graph}, + node_budget=15, + root=0, + ) + with pytest.raises(RuntimeError, match="before finding"): + mwcsb.get_solution(time_limit=1e-9) + + node_weighted = NodeWeightedSteinerProblem(graph, [[0, 15]], {v: 1 for v in graph}) + with pytest.raises(RuntimeError, match="before finding"): + node_weighted.get_solution(time_limit=1e-9) + + +@pytest.mark.parametrize("runner", ["core", "mwcsb"]) +def test_gurobi_tiny_positive_limit_is_status_safe(runner): + _require_gurobi() + graph = nx.cycle_graph(30) + nx.set_edge_attributes(graph, 1, "weight") + + if runner == "core": + problem = SteinerProblem(graph, [[0, 10, 20]], preprocess=False) + model, x, _y1, y2, z = mm.build_model_gurobi(problem, time_limit=1e-9) + result = mm.run_model_gurobi(model, problem, x, y2, z, return_status=True) + objective, edges = result[2], result[3] + else: + problem = BudgetedMaxWeightConnectedSubgraph( + graph, + {v: 1 for v in graph}, + {v: 1 for v in graph}, + node_budget=15, + root=0, + ) + model, _x, y1, _y2, _z, node_vars = mm.build_mwcsb_model_gurobi( + problem, time_limit=1e-9 + ) + result = mm.run_mwcsb_model_gurobi( + model, problem, y1, node_vars, return_status=True + ) + objective, edges = result[2], result[3] + + assert result[-1] in {"incomplete", "optimal"} + if result[-1] == "incomplete": + assert math.isinf(result[0]) + if not math.isfinite(objective): + assert edges == [] + + +def test_gurobi_mwcsb_infeasible_and_valid_unproven(monkeypatch): + _require_gurobi() + graph = nx.path_graph(3) + nx.set_edge_attributes(graph, 1, "weight") + + infeasible = BudgetedMaxWeightConnectedSubgraph( + graph, + {0: 5, 1: 1, 2: 1}, + {0: 2, 1: 1, 2: 1}, + node_budget=1, + root=0, + ) + model, _x, y1, _y2, _z, node_vars = mm.build_mwcsb_model_gurobi( + infeasible, time_limit=5 + ) + result = mm.run_mwcsb_model_gurobi( + model, infeasible, y1, node_vars, return_status=True + ) + assert result[-1] == "infeasible" + assert result[2] == -math.inf and result[3:5] == ([], []) + + feasible = BudgetedMaxWeightConnectedSubgraph( + graph, + {0: 5, 1: 1, 2: 10}, + {0: 0, 1: 1, 2: 1}, + node_budget=2, + root=0, + ) + model, _x, y1, _y2, _z, node_vars = mm.build_mwcsb_model_gurobi( + feasible, time_limit=5 + ) + monkeypatch.setattr( + mm, + "_gurobi_outcome", + lambda solved_model, _grb: ( + "incomplete", + solved_model.SolCount > 0, + ), + ) + result = mm.run_mwcsb_model_gurobi( + model, feasible, y1, node_vars, return_status=True + ) + assert result[-1] == "incomplete" + assert math.isinf(result[0]) and math.isfinite(result[2]) + assert set(result[4]) == {0, 1, 2} diff --git a/tests/test_steiner_problem.py b/tests/test_steiner_problem.py index 2665dc2..396b491 100644 --- a/tests/test_steiner_problem.py +++ b/tests/test_steiner_problem.py @@ -1546,10 +1546,11 @@ def test_solve_sap_highs_exhausted_time_limit_gap(): ctx = PrizeCollectingProblem(g, [[0]], prizes, penalty_cost=0)._build_pc_transform() view = DirectedSteinerProblem(ctx.sap_graph, ctx.root, ctx.terminals, weight=ctx.weight) - # With a dual-ascent upper bound, report an honest gap against it... + # An upper bound alone is not enough to certify a gap when no solve ran: + # there is no valid lower-bound result to pair it with. gap_ub, *_ = solve_sap_highs(view, time_limit=0.0, da_ub=10.0) - assert gap_ub == pytest.approx(1.0) # (10 - 0) / max(1, 10) - # ...without one, the relaxation gap would be spurious, so report inf. + assert gap_ub == float("inf") + # Without one the result is likewise explicitly unknown. gap_none, *_ = solve_sap_highs(view, time_limit=0.0) assert gap_none == float("inf")