Coding Patterns, TDD & Debugging

Algorithm Engineering — T06

Dr. Dominik Krupke

Coding Patterns, TDD & Debugging

Optimization code that isn’t tested is optimization code you can’t trust.

How much structure does the code deserve?

Coding Principles

Software engineering principles carry over but with different emphasis.

Fail fast on bad data (or take the blame for it)

Bad input should fail early. Bad output should fail loudly. Python’s duck typing can hide bad data. Frameworks like Pydantic can enforce a schema and raise loudly when it is violated.

class KnapsackInstance(BaseModel):
    # PositiveInt enforces that weights, values, and capacity are all positive integers.
    weights: list[PositiveInt] = Field(..., description="Weight in KG of each item.")
    values: list[PositiveInt] = Field(..., description="Value in Euro of each item.")
    capacity: PositiveInt = Field(..., description="Maximum weight in KG the knapsack can carry.")

    @model_validator(mode="after")
    def check_lengths(cls, v):
        if len(v.weights) != len(v.values):
            raise ValueError("Mismatch in number of weights and values.")
        return v

    def __str__(self)-> str:  # For visual sanity check in the logs.
        return f"KnapsackInstance(weights={len(self.weights)}|min={min(self.weights)}|max={max(self.weights)}, \
                 values={len(self.values)}|min={min(self.values)}|max={max(self.values)}, capacity={self.capacity})"

Validate (and document) at the boundaries, whenever data changes ownership. There are many opportunities for critical bugs here.

Narrow interfaces: Feed the solver only what it needs

Failure mode: A database migration requires a refactor of the solver.

In production system, the data to work on often lives in the system’s database and it is too natural to let the solver just fetch what it needs from there, creating two problems:

The solver now has a strong dependency on the database schema, making any change risky.
You have to dig through details to figure out what the solver depends on.

Keep the instance schema narrow, explicit, and self-contained. Have a separate extraction layer.

Inspectable — where did the wrong coefficient come from?
Decision-revealing — what facts actually influence the solver’s decisions?
Serializable — dump it to JSON, replay and benchmark.
Change-robust — the database engineer does not need to look into the solver details.

Split Data Concerns

Not one giant mutable object threaded through the whole pipeline. Separate types for separate concepts.

ProblemInstance: the data describing the problem to be solved.
SolveConfig: solver/technology specific parameters, limits, and logging settings.
Solution: the data describing the solution returned by the solver.
SolutionStats: solver/technology-specific additions to the solution.

Three design tests for whether the split is right:

Can I have a pool of solutions and decide later which one is the best?
Can I carry-over a solution to a changed instance? (potentially needing repair)
Does my instance or solution restrict the technology?

For sequential optimization, things may split up even further.

Separate logical concerns as complexity grows

Do not start with unnecessary architecture. Once the code grows, separate responsibilities clearly:

Data layer: load, clean, map, validate input data.
Model layer: decision variables, constraints, objective.
Solver layer: solver choice, parameters, limits, logging.
Solution layer: extract a domain-level solution object.
Analysis layer: statistics, plots, reports, explanations.

A growth path, not a starting checklist: a simple knapsack earns none of these five layers.

Make the mathematics readable

Optimization code should expose the formulation, not bury it inside implementation detail.

accumulated_weight = sum(x * w for x, w in zip(xs, weights))
model.add(accumulated_weight <= capacity)

accumulated_value = sum(x * v for x, v in zip(xs, values))
model.maximize(accumulated_value)

A reader should spot, at a glance:

the decision variables,
the main constraints,
the objective,
and any important derived expressions, named, not inlined.

Modules and Hierarchy

Mathematical modeling often exploits block structure for decomposition techniques, but the same structure is also great for modularizing the code.

A small set of fundamental decisions forms the backbone; every other component only refines it.

class TourComp: ...                      # the backbone
class CapacityComp(TourComp): ...        # refines it
class TimeWindowComp(TourComp): ...      # refines it
class CostComp(TourComp, TimeWindowComp): # scoring logic

Components couple through the backbone, not to each other:

readable — from structure to details.
testable — isolate the components under test.
maintainable — swap/configure components.

Not every problem has such a structure. Knowing the problem abstractions often pays off here in identifying the structure.

Consider the solver as an oracle, not a calculator

Recompute the objective yourself. A mismatch with the solver log is a bug you just found.
- You often round or simplify the objective in the model anyway, so you need your own value regardless.
Recompute auxiliary values too, never read them off the model variables.
- Solver-agnostic postprocessing can be reused across solver variants.
Check feasibility independently, since numerical tolerances can pass slightly infeasible solutions.

Ideally, extract only the solution’s “genome” from the solver and compute everything else yourself. That makes it easiest to swap technologies and work at different abstraction levels; your postprocessing can even become a full simulation using the solver as an oracle.

Start development with solver logging enabled

solver.parameters.log_search_progress = True
status = solver.solve(model)

Starting CP-SAT solver v9.10.4067
Parameters: max_time_in_seconds: 30 log_search_progress: true

Initial optimization model '':
#Variables: 450 (#bools: 276 #ints: 6 in objective)
#kLinearN: 94 (#terms: 1'392)

Check the parameters line and the variable/constraint counts against what you meant to build, before looking at anything else in the log.

Know the complexities and capabilities of your solver

How a constraint is written, not just what it expresses, decides how hard the solver has to work.

Hand-rolled	Solver built-in	Why it wins
boolean multiplication for an AND	linear AND encoding	keeps the model linear, not quadratic
lexicographic objective by hand	native lexicographic support	exploits internals you cannot reach
big-M constraints	indicator constraints	faster and numerically more stable

The reverse also happens: a built-in is generic by design. With exploitable domain knowledge, a hand-tailored formulation can beat what the solver would rederive. Default to the cheapest built-in; override it only where you genuinely know something the solver cannot.

Save inputs and outputs for inspection

def add_test_case(instance, config):
    solution = solve_knapsack(instance, config)
    (test_folder / "instance.json").write_text(instance.model_dump_json())
    (test_folder / "config.json").write_text(config.model_dump_json())
    (test_folder / "solution.json").write_text(solution.model_dump_json())

When suspicious behavior appears, a saved instance and solution let you:

reproduce the exact case, not a description of it,
validate the saved solution against the independent checker from a few slides ago,
compare old and new implementations on the same input,
promote the incident directly into a permanent regression test.

A bug report without the triggering instance is often difficult to investigate. A bug report with the exact instance and solution is usually actionable the same day.

Coding Patterns

Six patterns, mostly on the same running example: knapsack.

Pattern 1: Single Function for Simple Use Cases

def solve_knapsack(
    weights: list[int], values: list[int], capacity: int
) -> list[int]:
    model = cp_model.CpModel()
    n = len(weights)
    x = [model.new_bool_var(f"x_{i}") for i in range(n)]

    accumulated_weight = sum(weights[i] * x[i] for i in range(n))
    model.add(accumulated_weight <= capacity)
    accumulated_value = sum(values[i] * x[i] for i in range(n))
    model.maximize(accumulated_value)

    solver = cp_model.CpSolver()
    status = solver.solve(model)
    ok = status in (cp_model.OPTIMAL, cp_model.FEASIBLE)
    return [i for i in range(n) if solver.value(x[i])] if ok else []

Every line of code is a liability. Only use the more complex patterns and principles if they are needed to manage complexity, testability, or change. Otherwise, a single function is the right amount of structure.

Pattern 2: Proper Input and Output Schemas

class KnapsackInstance(BaseModel):
    weights: list[PositiveInt] = Field(..., description="Weight of each item, must be positive")
    values: list[PositiveInt] = Field(..., description="Value of each item, must be positive")
    capacity: PositiveInt = Field(..., description="Maximum weight the knapsack can carry, must be positive")

class KnapsackSolverConfig(BaseModel):
    time_limit: PositiveFloat = Field(10.0, description="Time limit for the solver in seconds, must be positive")
    opt_tol: NonNegativeFloat = Field(0.0, description="Optimality tolerance for the solver, must be non-negative")

type ItemIdx = int

class KnapsackSolution(BaseModel):
    selected_items: list[ItemIdx] = Field(..., description="Indices of items selected to pack in the knapsack")

The schema must already be fully prepared for optimization. No preprocessing inside the model, that would combine two hard jobs into one unmaintainable one.

Pattern 3: An Index for Efficient and Readable Queries

Certain queries may require indices that do not belong into the pure data model.

class NurseRosteringIndex:
    """Built once from the instance; each query is a dict lookup, not a scan."""
    def __init__(self, instance: NurseRosteringInstance): ...

    def nurses_with_skill(self, skill: Skill) -> list[NurseId]: ...
    def available_for(self, shift: ShiftId) -> list[NurseId]: ...
    def qualified_for(self, shift: ShiftId) -> list[NurseId]: ...
    def feasible_shifts(self, nurse: NurseId) -> list[ShiftId]: ...
    def shifts_on(self, day: Day) -> list[ShiftId]: ...
    def conflicting_shifts(self, shift: ShiftId) -> list[ShiftId]: ...  # overlapping

for shift in instance.shifts:                              # covering constraint
    qualified = index.nurses_with_skill(shift.required_skill)
    model.add_at_least_one(x[n, shift.id] for n in qualified)  # no scan over all nurses

These indices are generally very nice to unit test. Move as much of the modelling logic into them as possible.

Pattern 4: Class Instead of Function

class KnapsackSolver:
    def __init__(self, instance: KnapsackInstance, config: KnapsackSolverConfig):
        self.instance, self.config = instance, config
        self.model = cp_model.CpModel()
        self.x = [self.model.new_bool_var(f"x_{i}") for i in range(len(instance.weights))]
        self._build_model()
        self.solver = cp_model.CpSolver()

    def solve(self, time_limit=None) -> KnapsackSolution:
        ...  # set parameters, solve, extract

    def prohibit_combination(self, item_a: int, item_b: int):
        """Discovered after showing the user a solution: these two must not co-pack."""
        self.model.add(self.x[item_a] + self.x[item_b] <= 1)

solver = KnapsackSolver(instance, config)
solution = solver.solve()
solver.prohibit_combination(0, 1)      # a downstream check reveals a new rule
solution = solver.solve(time_limit=5)  # re-solve, no rebuild

A stateless function cannot express “add this constraint and re-solve.” A class holding both model and solver can, without rebuilding from scratch.

Pattern 5: Variable Containers

class _ItemSelectionVars:
    def __init__(self, instance, model, var_name="x"):
        self.instance = instance
        self.x = [model.new_bool_var(f"{var_name}_{i}")
                  for i in range(len(instance.weights))]

    def packs_item(self, i):  return self.x[i]
    def used_weight(self):    return sum(w * xi for w, xi in zip(self.instance.weights, self.x))
    def packed_value(self):   return sum(v * xi for v, xi in zip(self.instance.values, self.x))

self.model.add(self._item_vars.used_weight() <= self.instance.capacity)
self.model.maximize(self._item_vars.packed_value())

Potentially lazy generation: mint a variable in __getitem__ on first access instead of up front, when only a few of quadratically many candidates are ever used.

Pattern 6: Submodels

from cpsat_utils.piecewise import PiecewiseLinearFunction

f_costs = PiecewiseLinearFunction.from_points([(0, 0), (1000, 400), (1500, 1300)])
f_gain  = PiecewiseLinearFunction.from_points([(0, 0), (100, 800), (300, 2000)])

y_cost = f_costs.add_lower_bound(model, buy_1)     # y >= f(buy_1)
y_gain = f_gain.add_upper_bound(model, produce_1)  # y <= f(produce_1)
model.maximize(y_gain - y_cost)                    # internal vars stay hidden

Hiding the internal auxiliary variables behind add_lower_bound/add_upper_bound buys:

reusable — the same submodel plugs into knapsack, nurse rostering, anything with the right shape,
swappable — a faster internal encoding later, without touching a single caller,
testable in isolation — the real prize, a tiny model with nothing to do with the surrounding problem.

Test-Driven Optimization

Writing tests early clarifies requirements before they cost you.

Scenario Builder Pattern

def test_back_to_back_shifts_violate_min_rest():
    (
        NurseScenarioBuilder()
        .add_shift(day=0, start_hour=0, length=8)
        .add_shift(day=0, start_hour=8, length=8)  # starts exactly when the first ends
        .add_nurse("Alice", min_time_between_shifts=timedelta(hours=8))
        .assign("Alice", shifts=[0, 1])             # forced onto both, no rest between them
        .assert_infeasible()
    )

The defining trait: each call mutates internal state and returns self, so calls chain, and the terminal call, here assert_infeasible(), is the one hidden step that builds the CP-SAT model, wires in every module, solves, and checks the outcome.

Test submodules in isolation

def run_min_rest_test(assignments, expected_feasible):
    shifts = create_shifts(len(assignments), shift_length=8)
    nurse = create_nurse("Nurse A", min_time_between_shifts=timedelta(hours=8))
    instance = NurseRosteringInstance(nurses=[nurse], shifts=shifts)
    context = AssertModelFeasible() if expected_feasible else AssertModelInfeasible()
    with context as model:
        nurse_vars = NurseDecisionVars(nurse, shifts, model)
        MinTimeBetweenShifts().build(instance, model, [nurse_vars])
        for shift, assign in zip(shifts, assignments):
            if assign is not None:
                nurse_vars.fix(shift.uid, assign)

run_min_rest_test([True, False, True], expected_feasible=True)        # every shift fixed
run_min_rest_test([None, True, True, None], expected_feasible=False)  # only the colliding pair fixed

Feasibility: fix every assignment you can and push it right to the limit. Infeasibility: fix only the shifts that actually collide, leave everything else free.

A submodule tested with no surrounding problem

from cpsat_utils.piecewise import PiecewiseLinearFunction
from cpsat_utils.testing import assert_objective

def test_add_upper_bound_convex():
    """Concave-shaped function: a single convex part, no auxiliary variables."""
    f = PiecewiseLinearFunction([0, 10, 20], [0, 10, 5])
    model = cp_model.CpModel()
    x = model.new_int_var(0, 20, "x")
    y = f.add_upper_bound(model, x)
    model.maximize(y)
    solver = assert_objective(model, 10)  # push x to the peak
    assert solver.value(x) == 10

def test_add_upper_bound_non_convex():
    """Non-convex function: needs reified segment-selector constraints internally."""
    f = PiecewiseLinearFunction([0, 10, 20], [0, 10, 50])
    model = cp_model.CpModel()
    x = model.new_int_var(0, 20, "x")
    y = f.add_upper_bound(model, x)
    model.maximize(y)
    solver = assert_objective(model, 50)
    assert solver.value(x) == 20

No knapsack, no nurse rostering, no domain in sight: the same maximize-and-check-the-peak discipline from the previous slide, now applied to a function with no surrounding problem at all.

`cpsat-utils` test helpers

from cpsat_utils.testing import (
    AssertModelFeasible, AssertModelInfeasible, assert_objective,
)

with AssertModelFeasible() as model:
    ...   # raises RuntimeError if the model turns out infeasible

with AssertModelInfeasible() as model:
    ...   # raises RuntimeError if the model turns out feasible

assert_objective(model=model, expected=-1.0)  # feasible AND optimum matches

Prove the test bites before it passes

class NoBlockedShiftsModule(ShiftAssignmentModule):
    def build(self, instance, model, nurse_shift_vars):
        return 0  # TODO: nothing enforced yet

def test_no_blocked_shifts_infeasible():
    """A nurse forced onto a blocked shift must make the model infeasible."""
    with AssertModelInfeasible() as model:
        nurse_vars = NurseDecisionVars(nurse, shifts, model)
        NoBlockedShiftsModule().build(instance, model, [nurse_vars])
        nurse_vars.fix(shifts[0].uid, True)  # forcing the blocked assignment

Immune to symmetries, not to fixed seeds

new_solution = solve_knapsack(instance, config)
assert new_solution.objective <= solution.upper_bound
assert solution.objective <= new_solution.upper_bound
# Do not test for the selected items: the solver might return a
# different solution of the same quality.

Property-based testing

from hypothesis import given, strategies as st

@given(st.lists(st.integers()))
def test_sort_properties(xs):
    ys = sorted(xs)
    assert len(ys) == len(xs)                       # length preserved
    assert all(a <= b for a, b in zip(ys, ys[1:]))   # non-decreasing
    assert Counter(ys) == Counter(xs)                # same multiset

Debugging & Explaining Infeasibility

No error, and no correct solution either. Where is the bug?

Debugging a model

By hand

Turn every constraint off, then switch them back on one at a time, re-solving after each. The constraint that flips the model from fine to broken points at the faulty component.

Often the fault is not a true mistake in one component but the interference of two: for example, two individually-correct symmetry-breaking rules that are incompatible together. A one-by-one search can miss it, since each constraint looks fine alone.

By asking the solver why

Running example: 2-coloring

Give each node one of two colors so adjacent nodes differ. The triangle A–B–C is an odd cycle, so no 2-coloring exists: the model is UNSAT.

Blame everything, or pinpoint the conflict

“The set of all constraints cannot be satisfied.” Not useful.

“This small subset cannot be satisfied together.”

Minimal Unsatisfiable Subset (MUS)

A subset \(C' \subseteq C\) is a MUS of \(C\) if:

\(C'\) is unsatisfiable: \(\text{solve}(C') = \text{UNSAT}\)
\(C'\) is minimal: \(\forall c \in C',\ \text{solve}(C' \setminus \{c\}) = \text{SAT}\)

Same object, different vocabulary: MIP solvers call this an IIS (irreducible infeasible subset).

Counterfactual: what do I remove to fix it?

A Minimal Correction Subset (MCS) is minimal such that \(C \setminus C'\) is satisfiable, the complement of a maximal satisfiable subset.

Drop the highlighted constraint and the graph becomes 2-colorable

Wrap-up

Optimization code that isn’t tested is optimization code you can’t trust.

Four movements, one habit

Good optimization code: narrow interfaces, readable math, a solver-independent spec.
Coding patterns: the concrete shapes, functions, classes, containers, submodels.
Test-driven optimization: write the test, watch it fail, make it pass.
Debugging and explainability: when it breaks, isolate first; when it is UNSAT, ask for a small subset, not a verdict on everything.

References

Guns, Tias, and Dimos Tsouros. 2024. L03: Solving, Debugging and Explanation Techniques. In: An Open-Source Course on Constraint Solving, KU Leuven. https://github.com/tias/constraint-solving-course.

Krupke, Dominik. 2024a. Cpsat-Utils: Testing and Modeling Helpers for Google’s OR-Tools CP-SAT Solver. Https://pypi.org/project/cpsat-utils/.

Krupke, Dominik. 2024b. The CP-SAT Primer. https://github.com/d-krupke/cpsat-primer.

Coding Patterns, TDD & Debugging

Coding Patterns, TDD & Debugging

How much structure does the code deserve?

Coding Principles

Fail fast on bad data (or take the blame for it)

Narrow interfaces: Feed the solver only what it needs

Split Data Concerns

Separate logical concerns as complexity grows

Make the mathematics readable

Modules and Hierarchy

Consider the solver as an oracle, not a calculator

Start development with solver logging enabled

Know the complexities and capabilities of your solver

Save inputs and outputs for inspection

Coding Patterns

Pattern 1: Single Function for Simple Use Cases

Pattern 2: Proper Input and Output Schemas

Pattern 3: An Index for Efficient and Readable Queries

Pattern 4: Class Instead of Function

Pattern 5: Variable Containers

Pattern 6: Submodels

Test-Driven Optimization

Scenario Builder Pattern

Test submodules in isolation

A submodule tested with no surrounding problem

cpsat-utils test helpers

Prove the test bites before it passes

Immune to symmetries, not to fixed seeds

Property-based testing

Debugging & Explaining Infeasibility

Debugging a model

Running example: 2-coloring

Blame everything, or pinpoint the conflict

Minimal Unsatisfiable Subset (MUS)

Counterfactual: what do I remove to fix it?

Wrap-up

Four movements, one habit

Further reading

References

`cpsat-utils` test helpers