Solution Validation Testing

Solution Validation and Testing

You are an expert in validating optimization code: establishing, before anyone else asks, that reported solutions are feasible, that objectives are computed correctly, and that claimed improvements are real. This skill covers the five-layer validation stack — independent feasibility checkers, objective recomputation separate from the model, unit tests for constraint builders and operators, known-optimum regression tests, and exact-vs-heuristic cross-validation on small instances. Use the protocols below whenever optimization results will be trusted, published, or acted on. Hooker (1995, "Testing Heuristics: We Have It All Wrong") and Johnson (2002, "A Theoretician's Guide to the Experimental Analysis of Algorithms") both observed that most reported algorithm comparisons fail at exactly this layer: the code computing the numbers was never independently checked.

Initial Assessment

Establish these facts before writing any validation code. The answers determine which layers of the stack you need and how deep each must go.

• Map every place a solution is born or transformed. Solver extraction (.X values), construction heuristics, decoders, crossover/mutation operators, repair routines, local-search moves, file deserialization. Each is a validation boundary; a checker call belongs at every one of them.
• List the constraint families from the problem statement, not from the model code. The checker must be derived from the written specification. If the only description of a constraint is the model code itself, write the prose statement first — that document is what the checker implements.
• Pin down the objective convention exactly. Minimization or maximization; units; rounding rules (TSPLIB EUC_2D rounds each distance to the nearest integer — Reinelt 1991, "TSPLIB — A Traveling Salesman Problem Library"); whether soft-constraint penalties are part of the reported objective or reported separately.
• Inventory the ground-truth sources. Brute-forceable instance sizes, published optima (OR-Library — Beasley 1990; MIPLIB; TSPLIB; QAPLIB), a second solver, a trusted prior implementation. No ground truth at any size means metamorphic testing carries more weight (see Advanced Techniques).
• Record the solver tolerances in play. Gurobi defaults: IntFeasTol 1e-5, FeasibilityTol 1e-6, MIPGap 1e-4. A solution the solver calls optimal can violate exact integrality and exact feasibility by these amounts; the checker must handle that boundary deliberately.
• Separate hard constraints from soft ones. Hard violations make a solution invalid; soft violations are costs. The checker must report them differently, and the objective recomputation must include penalty terms with the documented weights.
• Check for shared code between pipelines. If the heuristic, the exact model, and the "checker" all call the same objective function, one bug makes all three agree. Shared code is the single most common cause of false validation confidence.
• Establish reproducibility of failures. Are instances generated with seeds? Is every run's (instance, seed, parameters) tuple recorded? A validation failure that cannot be reproduced cannot be fixed.
• Set the compute budget per test tier. Per-commit suites should run in seconds without a solver license; nightly suites can afford small exact solves. Decide the budget now so the suite is actually run.
• Determine what is already tested. Existing pytest suites, ad-hoc assertion scripts, manual spot checks. Fold them into the stack rather than duplicating them.
• Assess the cost of a wrong answer. A paper result, a production planning decision, or a thesis chapter each justify the full stack; a throwaway prototype may justify only layers 1-2. Scale effort to consequence, but never skip layer 1.

The Validation Stack

Five layers, ordered by how early they catch a bug and how much machinery they need. Build them in this order; each layer assumes the ones below it exist.

Formally, a checker is a function built only from the raw instance data $D$ and the written problem statement:

$$ \text{check}: (D, s) ;\mapsto; \big(, \text{feasible}(s) \in {\text{true}, \text{false}},;; V(s),;; \hat{f}(s) ,\big) $$

where $V(s)$ is the list of specific violations (never just a boolean) and $\hat{f}(s)$ is the independently recomputed objective. The independence rule: the only inputs the checker shares with the solver or heuristic are $D$ and the specification. No shared functions, no shared distance matrix builder, no shared penalty code. The checker is the specification in executable form, so it should be written in the most boring style available — plain loops, no cleverness, no performance tricks. Slow and obviously correct beats fast and probably correct here, because the checker is the court of last appeal.

Objective reconciliation uses a mixed absolute/relative tolerance:

$$ |f_{\text{solver}}(s) - \hat{f}(s)| ;\le; \varepsilon_{\text{abs}} + \varepsilon_{\text{rel}} \cdot |\hat{f}(s)| $$

with $\varepsilon_{\text{abs}} = \varepsilon_{\text{rel}} = 10^{-6}$ as a sane default for well-scaled data. Anything looser hides bugs; anything tighter generates false alarms from solver tolerances.

One asymmetry deserves emphasis because it drives the cross-validation protocol: a heuristic that beats a proven optimum has not found a better solution — it has found a bug. Either the heuristic's solution is infeasible (heuristic bug), or the exact model solves a different problem than the heuristic does (model bug). The independent checker is the arbiter: validate both solutions against the raw data and the side whose solution fails is the side that is wrong.

Decision tree

START: do you trust this optimization result?
|
+-- New project, no validation yet?
|     -> Build layer 1 (checker) and layer 2 (recomputation) BEFORE the first
|        real experiment. Add layer 3 alongside the first constraint builder
|        or operator. Add layers 4-5 as soon as a tiny instance solves.
|
+-- One solution looks suspicious?
|     -> Run layers 1+2 on it.
|        Infeasible        -> trace the violated family to the builder,
|                             decoder, or operator that owns it (layer 3).
|        Objective mismatch-> diff the objective term by term; check sense,
|                             rounding convention, penalty weights.
|
+-- Heuristic reports a better objective than the exact method?
|     -> Validate BOTH solutions with layers 1+2 against raw data.
|        Heuristic solution infeasible -> heuristic bug: missing constraint
|                                         check, incomplete repair.
|        Heuristic solution feasible   -> the exact model is wrong or solves
|                                         a different problem: audit each
|                                         constraint builder against the spec.
|
+-- Refactor, dependency upgrade, or parameter change planned?
|     -> Layer 4 regression suite green before AND after. Pin objective
|        values, never solution vectors (ties and symmetry break the latter).
|
+-- Results headed into a paper or a production decision?
      -> Full layer 5 sweep on small instances, layers 1+2 on every reported
         solution, and archive the validation table next to the result table.

Symptom-to-layer diagnosis

Layers 1 and 2 — Independent Feasibility Checking and Objective Recomputation

Design rules for checkers, in order of importance:

1. Read only raw instance data. Coordinates, demands, capacities — not the model's preprocessed arrays, not the heuristic's distance matrix. If preprocessing has a bug, a checker that consumes preprocessed data inherits it.
2. Collect every violation with a specific message. A boolean checker tells you that something is wrong; a violation list tells you what to fix and which component owns the bug. Do not stop at the first violation.
3. Recompute the objective even when infeasible. The objective of an infeasible solution is still diagnostic — a heuristic that is 2% better and slightly infeasible is a constraint-handling problem, not an objective problem.
4. Keep it deliberately simple. Plain Python loops over the specification. The checker for a problem with an $O(n)$ evaluation may be $O(n^2)$; it runs on demand, not in the inner loop.
5. Return a structured report, so harnesses (layer 5) and CI can consume it without parsing strings.

The end-to-end example below is a CVRP checker — a good template because the CVRP has several independent constraint families (vehicle count, visit-exactly-once, capacity) plus a distance objective with an explicit convention.

"""Independent feasibility checker and objective recomputation for the CVRP.

The checker reads only raw instance data (coordinates, demands, capacity).
It never imports the model or heuristic code, so a bug there cannot hide here.
"""
from dataclasses import dataclass, field
import math


@dataclass(frozen=True)
class CVRPInstance:
    """Raw CVRP data. Node 0 is the depot."""
    coords: tuple[tuple[float, float], ...]   # (x, y) per node, depot first
    demands: tuple[int, ...]                  # demand per node, demands[0] == 0
    capacity: int
    n_vehicles: int


@dataclass
class ValidationReport:
    """Outcome of one validation pass."""
    feasible: bool
    objective: float
    violations: list[str] = field(default_factory=list)

    def add(self, message: str) -> None:
        """Record a violation and mark the report infeasible."""
        self.feasible = False
        self.violations.append(message)


def euclid(a: tuple[float, float], b: tuple[float, float]) -> float:
    """Plain Euclidean distance — the documented objective convention."""
    return math.hypot(a[0] - b[0], a[1] - b[1])


def check_cvrp(instance: CVRPInstance, routes: list[list[int]]) -> ValidationReport:
    """Validate routes against raw data and recompute the objective.

    Each route is a customer sequence without the depot, e.g. [3, 1, 4].
    Checks: vehicle count, valid node ids, capacity, visit-exactly-once.
    """
    report = ValidationReport(feasible=True, objective=0.0)
    n = len(instance.coords)

    if len(routes) > instance.n_vehicles:
        report.add(f"{len(routes)} routes used, only "
                   f"{instance.n_vehicles} vehicles available")

    seen: dict[int, int] = {}
    for r_idx, route in enumerate(routes):
        load = 0
        for node in route:
            if not 1 <= node < n:
                report.add(f"route {r_idx}: invalid node id {node}")
                continue
            seen[node] = seen.get(node, 0) + 1
            load += instance.demands[node]
        if load > instance.capacity:
            report.add(f"route {r_idx}: load {load} exceeds "
                       f"capacity {instance.capacity}")
        path = [0, *route, 0]
        report.objective += sum(
            euclid(instance.coords[u], instance.coords[v])
            for u, v in zip(path, path[1:])
        )

    for customer in range(1, n):
        count = seen.get(customer, 0)
        if count != 1:
            report.add(f"customer {customer} visited {count} times, expected 1")
    return report


def reconcile_objective(model_obj: float, checker_obj: float,
                        abs_tol: float = 1e-6, rel_tol: float = 1e-6) -> bool:
    """True when the solver objective matches the recomputation within tolerance."""
    return abs(model_obj - checker_obj) <= abs_tol + rel_tol * abs(checker_obj)


# Tiny synthetic instance: depot at the origin, 4 customers.
inst = CVRPInstance(
    coords=((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (2.0, 0.0)),
    demands=(0, 3, 3, 3, 3),
    capacity=6,
    n_vehicles=2,
)
good = check_cvrp(inst, [[1, 2], [3, 4]])
bad = check_cvrp(inst, [[1, 2, 3], [4]])     # route 0 carries 9 > 6
print(good.feasible, round(good.objective, 4), bad.feasible, bad.violations[0])
# Expected: True 7.4142 False route 0: load 9 exceeds capacity 6

The second half of layers 1-2 is the boundary where solutions leave a MIP solver. Two rules: never read .X without checking Status and SolCount first, and never trust raw .X values of integer variables — they can sit IntFeasTol away from an integer. Round, measure the drift, and refuse to round when drift is large (large drift signals a numerically unstable model, not a rounding problem — see Klotz & Newman 2013, "Practical Guidelines for Solving Difficult Mixed Integer Linear Programs").

"""Extract a rounded solution from a solved Gurobi model, then re-verify it
from raw data. Never trust raw .X values of integer variables."""
import gurobipy as gp
from gurobipy import GRB


def extract_binary_solution(model: gp.Model, x: gp.tupledict) -> dict:
    """Round binary variables to {0,1} after checking solution availability."""
    if model.Status not in (GRB.OPTIMAL, GRB.TIME_LIMIT) or model.SolCount == 0:
        raise RuntimeError(f"no solution to extract (status {model.Status})")
    solution = {key: int(round(var.X)) for key, var in x.items()}
    drift = max(abs(var.X - round(var.X)) for var in x.values())
    if drift > 1e-4:  # far beyond default IntFeasTol=1e-5: investigate scaling
        raise RuntimeError(f"integrality drift {drift:.2e} too large to round")
    return solution


# Tiny demonstration: 0-1 knapsack, extract, then recheck from raw data.
values = [10, 13, 7]
weights = [4, 5, 3]
cap = 8
m = gp.Model("knapsack")
m.Params.OutputFlag = 0
y = m.addVars(3, vtype=GRB.BINARY, name="y")
m.addConstr(gp.quicksum(weights[i] * y[i] for i in range(3)) <= cap,
            name="capacity")
m.setObjective(gp.quicksum(values[i] * y[i] for i in range(3)), GRB.MAXIMIZE)
m.optimize()
sol = extract_binary_solution(m, y)
chosen = [i for i in range(3) if sol[i] == 1]
recomputed = sum(values[i] for i in chosen)        # independent of the model
weight_ok = sum(weights[i] for i in chosen) <= cap
print(chosen, recomputed, weight_ok, abs(m.ObjVal - recomputed) < 1e-6)
# Expected: [1, 2] 20 True True

Layer 3 — Unit Tests for Constraint Builders and Operators

Constraint builders (the add_<family>_constraints functions a well-structured model is made of — see milp-modeling-gurobi) are ordinary functions and deserve ordinary unit tests. The most effective strategy is probing: fix the variables to a hand-built solution with known feasibility status, then ask the solver whether the fixed model is feasible. This tests the builder against the problem statement, not against itself. Complement probing with structural assertions: constraint counts, constraint names, and objective coefficients read back from the model.

The running example is the generalized assignment problem (GAP): assign each job to exactly one machine at minimum cost, subject to machine capacities. Save this block as gap_model.py; the test and harness blocks import it — the tests must import the code under test, never re-implement it.

"""Generalized assignment model in gurobipy, built from named constraint
builders. Save as gap_model.py. Builders are separate functions so unit
tests can exercise each constraint family in isolation."""
import gurobipy as gp
from gurobipy import GRB
import numpy as np


def add_assignment_constraints(model: gp.Model, x: gp.tupledict,
                               n_machines: int, n_jobs: int) -> None:
    """Each job goes to exactly one machine."""
    for j in range(n_jobs):
        model.addConstr(
            gp.quicksum(x[i, j] for i in range(n_machines)) == 1,
            name=f"assign[{j}]",
        )


def add_capacity_constraints(model: gp.Model, x: gp.tupledict,
                             a: np.ndarray, b: np.ndarray) -> None:
    """Total resource use on each machine stays within its capacity."""
    n_machines, n_jobs = a.shape
    for i in range(n_machines):
        model.addConstr(
            gp.quicksum(a[i, j] * x[i, j] for j in range(n_jobs)) <= b[i],
            name=f"capacity[{i}]",
        )


def build_gap_model(c: np.ndarray, a: np.ndarray,
                    b: np.ndarray) -> tuple[gp.Model, gp.tupledict]:
    """Minimize assignment cost s.t. one-machine-per-job and capacities."""
    n_machines, n_jobs = c.shape
    model = gp.Model("gap")
    model.Params.OutputFlag = 0
    x = model.addVars(n_machines, n_jobs, vtype=GRB.BINARY, name="x")
    add_assignment_constraints(model, x, n_machines, n_jobs)
    add_capacity_constraints(model, x, a, b)
    model.setObjective(
        gp.quicksum(c[i, j] * x[i, j]
                    for i in range(n_machines) for j in range(n_jobs)),
        GRB.MINIMIZE,
    )
    return model, x


if __name__ == "__main__":
    # Tiny synthetic instance: 2 machines, 3 jobs.
    c = np.array([[4.0, 2.0, 5.0], [3.0, 6.0, 1.0]])
    a = np.array([[3.0, 2.0, 4.0], [2.0, 4.0, 2.0]])
    b = np.array([5.0, 5.0])
    model, x = build_gap_model(c, a, b)
    model.optimize()
    assignment = dict(sorted((j, i) for (i, j) in x if x[i, j].X > 0.5))
    print(model.Status == GRB.OPTIMAL, round(model.ObjVal, 1), assignment)
# Expected: True 6.0 {0: 1, 1: 0, 2: 1}

The probing tests. Note DualReductions = 0 so an infeasible fixed model reports INFEASIBLE rather than the ambiguous INF_OR_UNBD.

"""Unit tests for the GAP constraint builders. Save as test_gap_builders.py
and run with: pytest test_gap_builders.py -q

Strategy: probe each constraint family with hand-built assignments whose
feasibility is known from the problem statement; add structural assertions
on counts, names, and objective coefficients."""
import gurobipy as gp
from gurobipy import GRB
import numpy as np

from gap_model import build_gap_model


def fix_and_solve(model: gp.Model, x: gp.tupledict,
                  assignment: dict[int, int]) -> int:
    """Fix x to the given job->machine map and return the Gurobi status."""
    model.Params.DualReductions = 0   # report INFEASIBLE, not INF_OR_UNBD
    for (i, j), var in x.items():
        value = 1.0 if assignment.get(j) == i else 0.0
        var.LB = value
        var.UB = value
    model.optimize()
    return model.Status


def small_instance() -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """2 machines x 3 jobs; capacities chosen so some assignments overload."""
    c = np.array([[4.0, 2.0, 5.0], [3.0, 6.0, 1.0]])
    a = np.array([[3.0, 2.0, 4.0], [2.0, 4.0, 2.0]])
    b = np.array([5.0, 5.0])
    return c, a, b


def test_feasible_assignment_accepted() -> None:
    """A hand-checked feasible assignment must keep the model feasible."""
    model, x = build_gap_model(*small_instance())
    status = fix_and_solve(model, x, {0: 1, 1: 0, 2: 1})  # loads: m0=2, m1=4
    assert status == GRB.OPTIMAL


def test_overload_rejected() -> None:
    """Assigning all jobs to machine 0 overloads it (3+2+4 = 9 > 5)."""
    model, x = build_gap_model(*small_instance())
    status = fix_and_solve(model, x, {0: 0, 1: 0, 2: 0})
    assert status == GRB.INFEASIBLE


def test_constraint_counts_and_names() -> None:
    """One assignment row per job, one capacity row per machine, all named."""
    model, _ = build_gap_model(*small_instance())
    model.update()
    names = {constr.ConstrName for constr in model.getConstrs()}
    assert {"assign[0]", "assign[1]", "assign[2]"} <= names
    assert {"capacity[0]", "capacity[1]"} <= names
    assert model.NumConstrs == 5


def test_objective_coefficients_match_data() -> None:
    """Every x[i,j] must carry cost c[i,j] in the objective."""
    c, a, b = small_instance()
    model, x = build_gap_model(c, a, b)
    model.update()
    for (i, j), var in x.items():
        assert var.Obj == c[i, j]
# Expected: 4 passed when run under pytest with gap_model.py importable.

Metaheuristic components need a different test style: seeded random property tests. Instead of one hand-picked case, generate hundreds of random inputs from np.random.default_rng(seed) and assert an invariant that must hold for every one of them — the approach popularized by QuickCheck (Claessen & Hughes 2000, "QuickCheck: A Lightweight Tool for Random Testing of Haskell Programs"; the hypothesis library brings it to Python with automatic shrinking of failing cases). The two invariants below catch the two most damaging metaheuristic bugs in practice: operators that silently corrupt the representation, and delta evaluation that drifts away from full evaluation.

"""Property tests for metaheuristic components: a permutation crossover and
2-opt delta evaluation. Seeded random cases, fully reproducible."""
import numpy as np


def order_crossover(p1: np.ndarray, p2: np.ndarray,
                    rng: np.random.Generator) -> np.ndarray:
    """OX: copy a random slice of p1, fill the rest in p2's order."""
    n = p1.size
    lo, hi = np.sort(rng.choice(n + 1, size=2, replace=False))
    child = -np.ones(n, dtype=int)
    child[lo:hi] = p1[lo:hi]
    fill = p2[~np.isin(p2, child[lo:hi])]
    child[np.flatnonzero(child < 0)] = fill
    return child


def tour_length(tour: np.ndarray, dist: np.ndarray) -> float:
    """Full O(n) recomputation — the trusted reference."""
    return float(dist[tour, np.roll(tour, -1)].sum())


def two_opt_delta(tour: np.ndarray, dist: np.ndarray, i: int, k: int) -> float:
    """O(1) length change of reversing tour[i+1..k] (requires i < k)."""
    n = tour.size
    a, b = tour[i], tour[(i + 1) % n]
    c, d = tour[k], tour[(k + 1) % n]
    return float(dist[a, c] + dist[b, d] - dist[a, b] - dist[c, d])


def test_ox_always_yields_permutation() -> None:
    """Invariant: OX output is a permutation of 0..n-1 for ANY parents."""
    rng = np.random.default_rng(42)
    for _ in range(500):
        n = int(rng.integers(2, 30))
        p1, p2 = rng.permutation(n), rng.permutation(n)
        child = order_crossover(p1, p2, rng)
        assert np.array_equal(np.sort(child), np.arange(n))


def test_delta_matches_full_recomputation() -> None:
    """Invariant: delta == new length - old length, within float tolerance."""
    rng = np.random.default_rng(7)
    for _ in range(500):
        n = int(rng.integers(4, 40))
        pts = rng.random((n, 2))
        dist = np.linalg.norm(pts[:, None, :] - pts[None, :, :], axis=2)
        tour = rng.permutation(n)
        i, k = np.sort(rng.choice(n, size=2, replace=False))
        delta = two_opt_delta(tour, dist, int(i), int(k))
        new_tour = tour.copy()
        new_tour[i + 1:k + 1] = new_tour[i + 1:k + 1][::-1]
        full = tour_length(new_tour, dist) - tour_length(tour, dist)
        assert abs(delta - full) <= 1e-9 * (1.0 + abs(full)), (i, k, delta, full)


test_ox_always_yields_permutation()
test_delta_matches_full_recomputation()
print("operator property tests passed")
# Expected: operator property tests passed

The same pattern covers the other component types: a repair operator's output must always pass the layer-1 checker; a decoder must map every random-key vector to a valid phenotype; a neighborhood move applied and then undone must restore the original solution bit for bit. Write one invariant test per component, 200-500 random cases each, all seeded.

Layers 4 and 5 — Known-Optimum Regression and Exact-vs-Heuristic Cross-Validation

Layer 4 needs ground truth. Three sources, in order of preference:

1. Brute force on tiny instances. Enumerate the full solution space. The reference solver must be more trustworthy than the code it judges, so it is written in the most naive style possible — no pruning, no cleverness. Sizes: permutations up to $n \approx 9$ ($9! \approx 3.6 \times 10^5$), machine assignments up to $m^n \approx 10^6$, subsets up to $n \approx 20$.
2. Published optima for standard benchmarks (OR-Library, TSPLIB, QAPLIB, MIPLIB). These also validate your instance parser — but only after you confirm the objective convention matches (see Practical Challenges).
3. A previous trusted version of your own pipeline, pinned by commit hash.

Pin the values as data in the test file. A pinned value is a contract: a failing regression test means behavior changed, and the change must be explained before the pin is updated.

"""Brute-force GAP reference solver, seeded generator, and pinned regression
tests. Save as gap_reference.py; run the test with: pytest gap_reference.py -q

The brute force enumerates every machine choice per job: m^n combinations.
Keep n_jobs <= 8. It is deliberately naive — the reference must be more
trustworthy than the code it judges."""
import itertools

import numpy as np


def make_gap_instance(seed: int, n_machines: int,
                      n_jobs: int) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Seeded random GAP instance; capacities sized to keep most feasible."""
    rng = np.random.default_rng(seed)
    c = rng.integers(1, 20, size=(n_machines, n_jobs)).astype(float)
    a = rng.integers(1, 10, size=(n_machines, n_jobs)).astype(float)
    b = np.full(n_machines, np.ceil(1.3 * a.mean() * n_jobs / n_machines))
    return c, a, b


def brute_force_gap(c: np.ndarray, a: np.ndarray,
                    b: np.ndarray) -> tuple[float, tuple[int, ...]]:
    """Exact optimum by full enumeration: (cost, machine-per-job tuple)."""
    n_machines, n_jobs = c.shape
    best_cost, best_assign = float("inf"), ()
    for assign in itertools.product(range(n_machines), repeat=n_jobs):
        load = np.zeros(n_machines)
        for j, i in enumerate(assign):
            load[i] += a[i, j]
        if np.any(load > b + 1e-9):
            continue
        cost = float(sum(c[i, j] for j, i in enumerate(assign)))
        if cost < best_cost:
            best_cost, best_assign = cost, assign
    return best_cost, best_assign


# Pinned fixtures: (seed, n_machines, n_jobs) -> optimum, computed once by
# the brute force above and frozen. A red test means the generator or the
# reference changed behavior; explain the change before updating a pin.
KNOWN_OPTIMA: dict[tuple[int, int, int], float] = {
    (0, 3, 6): 33.0,
    (7, 3, 6): 45.0,
    (11, 3, 6): 38.0,
    (23, 3, 6): 18.0,
}


def test_known_optima_still_reproduced() -> None:
    """Regression: the (generator, reference) pair reproduces pinned optima."""
    for (seed, n_machines, n_jobs), pinned in KNOWN_OPTIMA.items():
        c, a, b = make_gap_instance(seed, n_machines, n_jobs)
        cost, _ = brute_force_gap(c, a, b)
        assert cost == pinned, f"seed {seed}: got {cost}, pinned {pinned}"


if __name__ == "__main__":
    test_known_optima_still_reproduced()
    best, arg = brute_force_gap(*make_gap_instance(0, 3, 6))
    print(best, arg)
# Expected: 33.0 (1, 1, 0, 0, 2, 0)

Layer 5 closes the loop: run the exact model and the heuristic over a bed of seeded small instances and assert five invariants per instance. Each invariant, when it fails, points at a specific component:

The complete harness:

"""Cross-validation harness: exact (Gurobi, MIPGap=0) versus a greedy +
first-improvement heuristic on a bed of seeded small GAP instances.

Invariants per instance:
  I1 exact solution passes the independent checker
  I2 exact objective reconciles with the recomputation
  I3 heuristic solution passes the independent checker
  I4 heuristic objective >= optimum - tol   (minimization)
  I5 Gurobi optimum == brute-force optimum
"""
import gurobipy as gp
from gurobipy import GRB
import numpy as np
import pandas as pd

from gap_model import build_gap_model
from gap_reference import brute_force_gap, make_gap_instance


def check_gap(assign: dict[int, int], c: np.ndarray, a: np.ndarray,
              b: np.ndarray) -> tuple[bool, float, list[str]]:
    """Independent GAP checker: each job assigned once, capacities respected."""
    n_machines, n_jobs = c.shape
    violations: list[str] = []
    if sorted(assign) != list(range(n_jobs)):
        violations.append("not every job assigned exactly once")
    load = np.zeros(n_machines)
    cost = 0.0
    for j, i in assign.items():
        load[i] += a[i, j]
        cost += c[i, j]
    for i in np.flatnonzero(load > b + 1e-9):
        violations.append(f"machine {i}: load {load[i]:.1f} > "
                          f"capacity {b[i]:.1f}")
    return (not violations), cost, violations


def solve_exact(c: np.ndarray, a: np.ndarray, b: np.ndarray,
                time_limit: float = 30.0) -> tuple[float, dict[int, int]]:
    """Proven optimum via the gurobipy model."""
    model, x = build_gap_model(c, a, b)
    model.Params.MIPGap = 0.0
    model.Params.TimeLimit = time_limit
    model.optimize()
    if model.Status != GRB.OPTIMAL:
        raise RuntimeError(f"exact solve failed: status {model.Status}")
    assign = dict(sorted((j, i) for (i, j) in x if x[i, j].X > 0.5))
    return model.ObjVal, assign


def greedy_heuristic(c: np.ndarray, a: np.ndarray,
                     b: np.ndarray) -> dict[int, int]:
    """Cheapest-feasible-machine construction + first-improvement moves."""
    n_machines, n_jobs = c.shape
    load = np.zeros(n_machines)
    assign: dict[int, int] = {}
    for j in np.argsort(a.min(axis=0))[::-1]:        # biggest jobs first
        feasible = [i for i in range(n_machines)
                    if load[i] + a[i, j] <= b[i] + 1e-9]
        if not feasible:
            feasible = list(range(n_machines))       # overload: I3 will flag it
        i = min(feasible, key=lambda m: c[m, j])
        assign[int(j)] = i
        load[i] += a[i, j]
    improved = True
    while improved:
        improved = False
        for j in range(n_jobs):
            i_old = assign[j]
            for i_new in range(n_machines):
                if (c[i_new, j] < c[i_old, j]
                        and load[i_new] + a[i_new, j] <= b[i_new] + 1e-9):
                    load[i_old] -= a[i_old, j]
                    load[i_new] += a[i_new, j]
                    assign[j] = i_new
                    improved = True
                    break
    return assign


def cross_validate(seeds: list[int], n_machines: int,
                   n_jobs: int) -> pd.DataFrame:
    """Assert invariants I1-I5 on every seeded instance; one row per instance."""
    rows = []
    for seed in seeds:
        c, a, b = make_gap_instance(seed, n_machines, n_jobs)
        bf_cost, _ = brute_force_gap(c, a, b)
        if bf_cost == float("inf"):
            continue                                  # infeasible instance
        exact_obj, exact_assign = solve_exact(c, a, b)
        ok_e, cost_e, viol_e = check_gap(exact_assign, c, a, b)
        assert ok_e, f"seed {seed}: exact solution fails checker: {viol_e}"
        assert abs(exact_obj - cost_e) <= 1e-6 * (1 + abs(cost_e)), \
            f"seed {seed}: objective mismatch {exact_obj} vs {cost_e}"
        assert abs(exact_obj - bf_cost) <= 1e-6 * (1 + abs(bf_cost)), \
            f"seed {seed}: Gurobi {exact_obj} != brute force {bf_cost}"
        heur_assign = greedy_heuristic(c, a, b)
        ok_h, cost_h, viol_h = check_gap(heur_assign, c, a, b)
        assert ok_h, f"seed {seed}: heuristic infeasible: {viol_h}"
        assert cost_h >= exact_obj - 1e-6, \
            f"seed {seed}: heuristic {cost_h} beats proven optimum {exact_obj}"
        rows.append({"seed": seed, "optimum": bf_cost, "heuristic": cost_h,
                     "gap_pct": 100.0 * (cost_h - bf_cost) / bf_cost})
    return pd.DataFrame(rows)


if __name__ == "__main__":
    table = cross_validate(seeds=list(range(20)), n_machines=3, n_jobs=7)
    print(table[["gap_pct"]].describe().loc[["mean", "max"]].round(1))
    print(f"validated {len(table)} instances, all invariants hold")
# Expected: gap_pct mean 8.2 and max 53.5 on these seeds, then
# "validated 20 instances, all invariants hold".

Two reading notes on the harness. First, I4 is an invariant, not a quality target — it must hold always; the heuristic's average gap (here 8.2%) is a benchmark observation that belongs in a monitored experiment table, never in an assert. Second, when I4 fires, resist the urge to "fix" the heuristic: run the checker on both solutions first. In practice the exact model is the guilty party about as often as the heuristic.

Advanced Techniques

Metamorphic testing of whole pipelines

When no ground truth exists at realistic sizes, test relations between runs instead of absolute values (Chen, Cheung & Yiu 1998, "Metamorphic Testing: A New Approach for Generating Next Test Cases"; Segura et al. 2016, "A Survey on Metamorphic Testing"). Useful relations for optimization pipelines: scaling all costs by $k > 0$ scales a minimum by exactly $k$; permuting entity labels leaves the optimum unchanged; relaxing a constraint (more capacity, longer deadline) never worsens a minimum; adding a dominated option (a machine worse in every cost and resource) leaves the optimum unchanged. These relations hold at any instance size, so they also validate the pipeline where brute force cannot reach.

"""Metamorphic tests: transform the instance, predict the optimum's response.
No ground truth needed — the relation between two runs IS the oracle."""
import numpy as np

from gap_reference import brute_force_gap, make_gap_instance


def optimum(c: np.ndarray, a: np.ndarray, b: np.ndarray) -> float:
    """Reference optimum for small instances."""
    cost, _ = brute_force_gap(c, a, b)
    return cost


def test_cost_scaling() -> None:
    """Scaling all costs by k > 0 must scale the optimum by exactly k."""
    c, a, b = make_gap_instance(seed=5, n_machines=3, n_jobs=6)
    base = optimum(c, a, b)
    assert abs(optimum(3.0 * c, a, b) - 3.0 * base) <= 1e-9 * (1.0 + base)


def test_machine_relabeling() -> None:
    """Permuting machine indices (rows of c, a, b) leaves the optimum fixed."""
    c, a, b = make_gap_instance(seed=5, n_machines=3, n_jobs=6)
    perm = np.array([2, 0, 1])
    assert optimum(c[perm], a[perm], b[perm]) == optimum(c, a, b)


def test_capacity_relaxation_never_hurts() -> None:
    """Adding capacity can only keep or improve a minimization optimum."""
    c, a, b = make_gap_instance(seed=5, n_machines=3, n_jobs=6)
    assert optimum(c, a, b + 5.0) <= optimum(c, a, b) + 1e-9


test_cost_scaling()
test_machine_relabeling()
test_capacity_relaxation_never_hurts()
print("metamorphic tests passed")
# Expected: metamorphic tests passed

To apply the same relations at large scale, replace brute_force_gap with the heuristic plus a fixed seed: relabeling invariance then additionally certifies that the heuristic has no hidden index-order dependence.

Fault injection: testing the tests

A checker that has never reported a violation is itself untested. After writing the validation stack, deliberately break the system one fault at a time and confirm the stack goes red: drop one constraint family from the model, flip an inequality direction, perturb a single objective coefficient by 1%, make the repair operator skip its last item, off-by-one a loop bound in the decoder. Every injected fault that survives all five layers marks a real gap in the stack — usually a constraint family the checker forgot or a tolerance set too loose. This is manual mutation testing; mutmut and cosmic-ray automate the mutation step for Python. Run a fault-injection pass when the checker is first written and after every change to the checker itself.

Tolerance engineering at the solver boundary

Solver tolerances are the main source of false alarms — and of falsely silent checkers. Three rules. First, round-then-verify: convert integer variables to exact integers (asserting drift below 1e-4) before checking, then perform all feasibility arithmetic on the rounded values; with integral data the capacity check then becomes exact integer arithmetic with no epsilon at all. Second, mind big-M leakage: with FeasibilityTol $= 10^{-6}$ and $M = 10^6$, a "satisfied" constraint $y \le M z$ admits $y = 1$ while $z = 0$ rounds to zero — the checker must test the logical condition ("no flow unless open"), not the linear surrogate. Third, scale the data: keep coefficient magnitudes within a range of about $10^6$ between smallest and largest, or tighten IntFeasTol and FeasibilityTol and accept the slowdown (Klotz & Newman 2013). When a checker disagrees with the solver by amounts near the tolerance, the model is numerically fragile; fix scaling rather than loosening the checker.

Differential testing across formulations and solvers

Two independent implementations of the same problem must agree on the optimal value — never require agreement on the solution, which is generally non-unique. Useful pairs: weak vs strong MIP formulations of the same problem; a MIP and a CP-SAT model; gurobipy and HiGHS on identical LP files; a compact and a pattern-based formulation (cutting stock). The disagreement protocol mirrors I4: validate both incumbents with the layer-1 checker against raw data; the side whose incumbent fails is wrong; if both pass, compare the proven bounds — the side whose bound excludes the other's feasible incumbent has the broken formulation. Differential testing earns its cost when a model is refactored: keep the old formulation alive in the test suite for one release as the comparison partner,

…(truncated)