Callbacks and User-Defined Cuts

discopt’s Branch and Bound solver supports three callback types that let you interact with the search process: node callbacks for monitoring progress, lazy constraint callbacks for adding cutting planes on the fly, and incumbent callbacks for filtering candidate solutions.

Callbacks are a standard mechanism in mixed-integer solvers [Quesada and Grossmann, 1992] for implementing problem-specific logic without modifying the solver itself. Lazy constraints, in particular, are essential for problems where the full constraint set is exponentially large (such as subtour elimination in the TSP) and constraints are only generated as needed [Sawaya and Grossmann, 2005].

Note

Callbacks never affect the correctness of the solver’s internal math. The solver catches and logs any exceptions raised by user callbacks, so a buggy callback will not crash the solve.

import discopt
import discopt.modeling as dm
import numpy as np
from discopt.callbacks import CallbackContext, CutResult

Node Callback: Logging B&B Progress

The simplest callback type is the node callback, which is called after each batch of B&B nodes is processed. It therefore fires only on a model that actually reaches the Python spatial-B&B loop: a pure MILP is routed to HiGHS and a convex MINLP to OA, and neither opens a B&B node, so a node callback attached to one is never invoked. The example below is a nonconvex MINLP for that reason, and the cell asserts the callback fired rather than assuming it. It receives a CallbackContext with the current state of the search (node count, incumbent objective, best bound, gap, and elapsed time). This is useful for custom logging, progress bars, or early termination logic.

# A nonconvex MINLP -- one that genuinely reaches the spatial-B&B loop.
# (A pure MILP like `min c'x` over binaries is routed to HiGHS and closes with
# node_count == 0, so a node callback attached to it is never called at all.)
np.random.seed(3)
n = 6
m = discopt.Model("progress_demo")
x = m.continuous("x", shape=(n,), lb=-2, ub=2)
y = m.binary("y", shape=(n,))
Q = np.random.randn(n, n)
Q = (Q + Q.T) / 2  # symmetric and indefinite -> nonconvex

m.minimize(
    dm.sum(
        lambda i: dm.sum(lambda j: float(Q[i, j]) * x[i] * x[j], over=range(n)),
        over=range(n),
    )
    + dm.sum(y)
)
m.subject_to(dm.sum(x) >= 1, name="total")
for i in range(n):
    m.subject_to(x[i] <= 2 * y[i], name=f"link_{i}")

# Define a node callback that logs progress
log_entries = []


def progress_logger(ctx: CallbackContext, model: discopt.Model) -> None:
    entry = {
        "nodes": ctx.node_count,
        "incumbent": ctx.incumbent_obj,
        "bound": ctx.best_bound,
        "gap": ctx.gap,
        "time": f"{ctx.elapsed_time:.3f}s",
    }
    log_entries.append(entry)


result = m.solve(node_callback=progress_logger, time_limit=30)
print(f"Final: status={result.status}, obj={result.objective:.4f}, "
      f"nodes={result.node_count}")
print(f"Node callback fired {len(log_entries)} time(s). First few:")
for entry in log_entries[:4]:
    print(f"  {entry}")

# Probe-fired check: without this, a callback that is never invoked reads as a pass.
assert log_entries, "the node callback was never invoked"
assert result.node_count > 0, "this model must reach the spatial-B&B loop"

Final: status=optimal, obj=-54.4997, nodes=79
Node callback fired 29 time(s). First few:
  {'nodes': 3, 'incumbent': -54.49969349863937, 'bound': -70.17695071295672, 'gap': 0.28765771342748475, 'time': '0.338s'}
  {'nodes': 7, 'incumbent': -54.49969349863937, 'bound': -69.63598320993978, 'gap': 0.27773164837483544, 'time': '0.422s'}
  {'nodes': 13, 'incumbent': -54.499694413900656, 'bound': -68.66919037710602, 'gap': 0.2599922094167064, 'time': '0.466s'}
  {'nodes': 21, 'incumbent': -54.499694413900656, 'bound': -67.68342652479917, 'gap': 0.24190469786443214, 'time': '0.494s'}

Lazy Constraints: Subtour Elimination for TSP

Lazy constraints are the most powerful callback type. They allow you to add linear cutting planes during the solve, which is essential for problems with an exponential number of constraints.

The classic example is the Travelling Salesman Problem (TSP), where subtour elimination constraints (SECs) prevent solutions that form disconnected cycles. Rather than adding all \(O(2^n)\) SECs upfront, we add them lazily: whenever the solver finds an integer-feasible solution that contains a subtour, we add a cut that eliminates it [Westerlund and Pettersson, 1995].

The callback receives a CallbackContext and returns a list of CutResult objects. Each CutResult specifies the cut as a list of (variable, coefficient) pairs, a sense ("<=", ">=", or "=="), and a right-hand side value.

Note

This section previously demonstrated subtour elimination on a 4-city TSP, where a subtour is structurally impossible under the degree-2 constraints — so the callback returned [] every time, the committed output read optimal, cost 80.0, and the feature being demonstrated was never exercised. The cell below now counts its invocations and the cuts it generates, so the no-op is visible, and it also runs a 6-city instance that does produce a subtour.

Exercising it that way found #1365: a callback that returned any CutResult made the solve come back status="unknown" with objective=None, because the node whose callback produced a cut was marked infeasible and its subtree fathomed. Returning [] was unaffected, which is why a 4-city demonstration could never have caught it. Fixed; the 6-city cell below asserts the answer.

# Small 4-city TSP with lazy subtour elimination
# Cities: 0, 1, 2, 3
# Decision variables: x[i,j] = 1 if edge (i,j) is in the tour
n = 4
# Distance matrix (symmetric)
dist = np.array(
    [
        [0, 10, 15, 20],
        [10, 0, 35, 25],
        [15, 35, 0, 30],
        [20, 25, 30, 0],
    ],
    dtype=float,
)

m = discopt.Model("tsp_4city")

# Binary variables for edges (upper triangle only, i < j)
edges = {}
edge_vars = []
for i in range(n):
    for j in range(i + 1, n):
        var = m.binary(f"x_{i}_{j}")
        edges[(i, j)] = var
        edge_vars.append((i, j, var))

# Minimize total distance
m.minimize(sum(dist[i, j] * var for i, j, var in edge_vars))

# Degree constraints: each city has exactly 2 edges
for city in range(n):
    incident = []
    for i, j, var in edge_vars:
        if i == city or j == city:
            incident.append(var)
    m.subject_to(sum(incident) == 2, name=f"degree_{city}")


def find_subtours(sol_flat, edge_list, n_cities):
    """Find connected components (subtours) in the solution."""
    # Build adjacency from solution
    adj = {i: [] for i in range(n_cities)}
    for idx, (i, j, _) in enumerate(edge_list):
        if sol_flat[idx] > 0.5:
            adj[i].append(j)
            adj[j].append(i)

    visited = set()
    components = []
    for start in range(n_cities):
        if start in visited:
            continue
        component = set()
        stack = [start]
        while stack:
            node = stack.pop()
            if node in component:
                continue
            component.add(node)
            for nb in adj[node]:
                if nb not in component:
                    stack.append(nb)
        visited |= component
        components.append(component)
    return components


def subtour_callback(ctx: CallbackContext, model: discopt.Model):
    """Add subtour elimination cuts for any disconnected component."""
    components = find_subtours(ctx.x_relaxation, edge_vars, n)

    if len(components) <= 1:
        return []  # Single tour, no subtours

    cuts = []
    for comp in components:
        if len(comp) == n:
            continue  # Full tour
        # SEC: sum of edges within the subtour <= |S| - 1
        terms = []
        for i, j, var in edge_vars:
            if i in comp and j in comp:
                terms.append((var, 1.0))
        if terms:
            cuts.append(
                CutResult(
                    terms=terms,
                    sense="<=",
                    rhs=float(len(comp) - 1),
                )
            )
    return cuts


# Count what the callback actually does, so a no-op cannot read as a success.
cb_stats = {"calls": 0, "cuts": 0}
_raw_callback = subtour_callback


def subtour_callback(ctx: CallbackContext, model: discopt.Model):
    cb_stats["calls"] += 1
    cuts = _raw_callback(ctx, model)
    cb_stats["cuts"] += len(cuts)
    return cuts


result = m.solve(lazy_constraints=subtour_callback, time_limit=60)
print(f"Status: {result.status}")
print(f"Optimal tour cost: {result.objective}")
print(f"Callback invoked {cb_stats['calls']} time(s); cuts generated: {cb_stats['cuts']}")
if result.x is not None:
    print("Edges in tour:")
    for i, j, var in edge_vars:
        val = result.x[f"x_{i}_{j}"]
        if np.all(val > 0.5):
            print(f"  {i} -- {j}  (dist={dist[i, j]:.0f})")

# On 4 cities every degree-2 solution is a single Hamiltonian cycle, so NO subtour
# can arise and no cut is ever needed. That is why this instance returns a correct
# answer: it does not exercise lazy constraints at all.
assert cb_stats["calls"] > 0, "the lazy-constraint callback was never invoked"
assert cb_stats["cuts"] == 0, (
    "a 4-city degree-2 TSP cannot have a subtour; a cut here would be a surprise"
)

Status: optimal
Optimal tour cost: 80.0
Callback invoked 2 time(s); cuts generated: 0
Edges in tour:
  0 -- 1  (dist=10)
  0 -- 2  (dist=15)
  1 -- 3  (dist=25)
  2 -- 3  (dist=30)

Where the cuts actually fire — and where #1365 shows

Six cities arranged as two tight triangles make the subtour the cheap answer, so the callback has real work to do. Without it the solver returns two disjoint 3-cycles, which is not a tour; with it, six subtour-elimination cuts are generated and the solve returns the real tour.

A Hamiltonian cycle on this instance must cross between the triangles at least twice, so the optimum is 2 x 5 + 4 x 1 = 14, against the 6.0 the uncut model reports for a pair of 3-cycles. The cell asserts both.

n6 = 6
dist6 = np.full((n6, n6), 5.0)
np.fill_diagonal(dist6, 0.0)
for a, b in [(0, 1), (1, 2), (0, 2), (3, 4), (4, 5), (3, 5)]:
    dist6[a, b] = dist6[b, a] = 1.0

m6 = discopt.Model("tsp6")
ev6 = []
for i in range(n6):
    for j in range(i + 1, n6):
        ev6.append((i, j, m6.binary(f"x_{i}_{j}")))
m6.minimize(sum(dist6[i, j] * v for i, j, v in ev6))
for c in range(n6):
    m6.subject_to(sum(v for i, j, v in ev6 if c in (i, j)) == 2, name=f"degree6_{c}")


def edges_of(xdict):
    return [(i, j) for i, j, _ in ev6 if xdict[f"x_{i}_{j}"] > 0.5]


plain = m6.solve(time_limit=60)
print(f"no callback: obj={plain.objective}, edges={edges_of(plain.x)}")
print(f"  components: {[sorted(c) for c in find_subtours(np.array([plain.x[f'x_{i}_{j}'] for i, j, _ in ev6]).ravel(), ev6, n6)]}")

stats6 = {"calls": 0, "cuts": 0}


def sec6(ctx: CallbackContext, model: discopt.Model):
    stats6["calls"] += 1
    comps = find_subtours(ctx.x_relaxation, ev6, n6)
    if len(comps) <= 1:
        return []
    out = []
    for comp in comps:
        if len(comp) == n6:
            continue
        terms = [(v, 1.0) for i, j, v in ev6 if i in comp and j in comp]
        if terms:
            out.append(CutResult(terms=terms, sense="<=", rhs=float(len(comp) - 1)))
    stats6["cuts"] += len(out)
    return out


cut = m6.solve(lazy_constraints=sec6, time_limit=60)
print(f"with SECs:   status={cut.status}, obj={cut.objective}, "
      f"calls={stats6['calls']}, cuts={stats6['cuts']}")
assert stats6["cuts"] > 0, "this instance must generate subtour cuts"

# Hand-derived: any Hamiltonian cycle crosses between the triangles twice, so the
# optimum is 2*5 + 4*1 = 14. The uncut model's 6.0 is two disjoint 3-cycles.
assert plain.objective == 6.0, plain.objective
assert cut.status == "optimal" and abs(cut.objective - 14.0) < 1e-6, (
    f"status={cut.status} obj={cut.objective}"
)
tour_edges = np.array([cut.x[f"x_{i}_{j}"] for i, j, _ in ev6]).ravel()
assert len(find_subtours(tour_edges, ev6, n6)) == 1, (
    "the cut solve must return a single tour"
)

no callback: obj=6.0, edges=[(0, 1), (0, 2), (1, 2), (3, 4), (3, 5), (4, 5)]
  components: [[0, 1, 2], [3, 4, 5]]
with SECs:   status=optimal, obj=14.0, calls=5, cuts=6

Incumbent Callback: Filtering Solutions

The incumbent callback is called whenever the solver finds a new best integer-feasible solution. You can inspect the solution and return False to reject it. This is useful for enforcing problem-specific feasibility conditions that are difficult to express as algebraic constraints.

For example, you might reject solutions that violate a complex business rule, or solutions that fail an external simulation check.

Note

ctx.incumbent_obj is the previous incumbent, reported in the solver’s internal minimization sense. On the maximization below its optimum of 16.0 therefore shows up in the callback as -16.0. The solution dict passed as the third argument is in the model’s own variables and needs no such correction.

# Example: reject solutions where fewer than 2 items are selected
m = discopt.Model("filter_demo")
x = m.binary("x", shape=(4,))
profits = np.array([10.0, 6.0, 4.0, 2.0])
m.maximize(sum(profits[i] * x[i] for i in range(4)))
m.subject_to(sum(3 * x[i] for i in range(4)) <= 7, name="capacity")


accepted = []


def require_min_items(ctx, model, solution):
    """Reject solutions with fewer than 2 selected items."""
    x_val = solution["x"]
    n_selected = int(np.sum(x_val > 0.5))
    accepted.append(n_selected)
    if n_selected < 2:
        print(f"  Rejected: only {n_selected} item(s) selected")
        return False
    print(f"  Accepted: {n_selected} items selected, obj={ctx.incumbent_obj}")
    return True


result = m.solve(incumbent_callback=require_min_items, time_limit=30)
print(f"\nFinal: status={result.status}, obj={result.objective}")
if result.x is not None:
    print(f"Selected items: {np.where(result.x['x'] > 0.5)[0]}")
assert accepted, "the incumbent callback was never invoked"

# The rejecting branch, which the accept-everything rule above never reaches.
# Capacity 7 at cost 3 each admits at most 2 items, so "at least 3" rejects every
# candidate -- and the solve ends with no incumbent, which is the callback working.
rejected = []


def require_three_items(ctx, model, solution):
    n_selected = int(np.sum(solution["x"] > 0.5))
    rejected.append(n_selected)
    return n_selected >= 3


strict = m.solve(incumbent_callback=require_three_items, time_limit=30)
print(f"\nRejecting rule: {len(rejected)} candidate(s) seen, sizes {rejected}")
print(f"Result with every candidate rejected: status={strict.status}, obj={strict.objective}")
assert rejected, "the rejecting callback was never invoked"
assert strict.objective is None, "no solution should survive an always-reject rule"

  Accepted: 2 items selected, obj=None
  Accepted: 2 items selected, obj=-16.0

Final: status=optimal, obj=16.0
Selected items: [0 1]

Rejecting rule: 8 candidate(s) seen, sizes [2, 2, 2, 2, 2, 2, 2, 2]

Result with every candidate rejected: status=unknown, obj=None

API Summary

The callback API consists of three components in discopt.callbacks:

CallbackContext is a dataclass passed to all callbacks containing the current B&B state: node_count, incumbent_obj (or None), best_bound, gap (or None), elapsed_time, x_relaxation (the current node’s solution as a flat numpy array), and node_bound.

CutResult specifies a linear cut via terms (a list of (Variable, coefficient) tuples), sense ("<=", ">=", or "=="), and rhs (the right-hand side value). Cuts are converted to dense constraint vectors internally and added to the solver’s cut pool.

The three callback types are passed as keyword arguments to Model.solve():

• node_callback(ctx, model) -> None: called after each batch of nodes.

• lazy_constraints(ctx, model) -> list[CutResult]: called at integer-feasible nodes, returns cuts to add.

• incumbent_callback(ctx, model, solution) -> bool: called when a new incumbent is found, returns False to reject.

All callbacks are optional and can be combined freely [Smith and Pantelides, 1999].