User guide · 2.0 development

Callbacks

SHOTpy exposes SHOT’s callbacks through solver.registerCallback(locations, fn). A callback is registered for one or several locations in the solution process, given as SHOTpy.CallbackLocation values combined with | or as a list:

Location

Context class

Called

Actions

InteriorPointSearch

InteriorPointSearchContext

After ESH checks existing points or runs its interior-point search, even if no point was found

setInteriorPoints(points)

DualBoundUpdate

DualBoundUpdateContext

When SHOT can accept a dual bound; before and after a direct MIQCQP solve

setDualBound(value)

PrimalCandidateSearch

PrimalCandidateSearchContext

Before the first dual solve, when SHOT collects primal candidates, and during finalization

addPrimalSolution(point)

PrimalCandidateCheck

PrimalCandidateCheckContext

Before a primal candidate is checked

rejectCandidate()

NewPrimalSolution

NewPrimalSolutionContext

After a primal solution has been stored

–

BoundsUpdated

BoundsUpdatedContext

After a primal or dual bound changes

–

TerminationCheck

TerminationCheckContext

When SHOT checks whether to terminate

–

HyperplaneSelection

HyperplaneSelectionContext

When SHOT selects the hyperplanes to add

addHyperplane(hyperplane)

BoundTightening

BoundTighteningContext

In setProblem(), after the bound tightening of the original and of the reformulated problem

setVariableBounds(index, lower, upper)

The callback is called as fn(ctx), where ctx is of the context class for the location. All contexts have the state of the solver, e.g., iterationNumber, dualBound, primalBound, relativeGap, elapsedTime, solutionStatistics, originalProblem and reformulatedProblem, and ctx.terminate() stops SHOT at its next termination check. Values that are not available yet, e.g., the dual bound before the first dual problem has been solved, are infinite.

The context can only be used while the callback runs: copy what you want to keep. If a callback raises an exception, SHOT terminates and solveProblem() raises it.

import SHOTpy

# ── Build the ex1223b MINLP ────────────────────────────────────────────────
solver = SHOTpy.Solver()

problem = SHOTpy.Problem(solver)
problem.name = "ex1223b_callbacks"

x1 = problem.addVariable("x1", SHOTpy.VariableType.Real, 0.0, 10.0)
x2 = problem.addVariable("x2", SHOTpy.VariableType.Real, 0.0, 10.0)
x3 = problem.addVariable("x3", SHOTpy.VariableType.Real, 0.0, 10.0)
b4 = problem.addVariable("b4", SHOTpy.VariableType.Binary)
b5 = problem.addVariable("b5", SHOTpy.VariableType.Binary)
b6 = problem.addVariable("b6", SHOTpy.VariableType.Binary)
b7 = problem.addVariable("b7", SHOTpy.VariableType.Binary)

problem.setObjective((b4 - 1.0)**2 + (b5 - 2.0)**2 + (b6 - 1.0)**2 - SHOTpy.log(1.0 + b7)
                 + (x1 - 1.0)**2 + (x2 - 2.0)**2 + (x3 - 3.0)**2, SHOTpy.ObjectiveDirection.Minimize)

problem.addConstraint(x1 + x2 + x3 + b4 + b5 + b6 <= 5, "e1")
problem.addConstraint(b6**2 + x1**2 + x2**2 + x3**2 <= 5.5, "e2")
problem.addConstraint(x1 + b4 <= 1.2, "e3")
problem.addConstraint(x2 + b5 <= 1.8, "e4")
problem.addConstraint(x3 + b6 <= 2.5, "e5")
problem.addConstraint(x1 + b7 <= 1.2, "e6")
problem.addConstraint(b5**2 + x2**2 <= 1.64, "e7")
problem.addConstraint(b6**2 + x3**2 <= 4.25, "e8")
problem.addConstraint(b5**2 + x3**2 <= 4.64, "e9")
problem.finalize()

# ── Example 1: track every new primal solution ────────────────────────────
primal_log = []

def on_new_primal(ctx):
    primal_log.append({
        "iter":      ctx.iterationNumber,
        "obj":       ctx.objectiveValue,
        "dual":      ctx.dualBound,
        "gap (%)":   round(ctx.relativeGap * 100, 4),
    })

solver.registerCallback(SHOTpy.CallbackLocation.NewPrimalSolution, on_new_primal)

# ── Example 2: stop after the relative gap is tight enough ───────────────
def early_stop(ctx):
    # Terminate once the gap drops below 1 %
    if ctx.relativeGap < 0.01:
        ctx.terminate()

solver.registerCallback(SHOTpy.CallbackLocation.TerminationCheck, early_stop)

# ── Suppress console output and solve ────────────────────────────────────
solver.updateSetting("Output.Console.LogLevel", 6)  # 6 = Off
assert solver.setProblem(problem)
assert solver.solveProblem()
# ── Results ───────────────────────────────────────────────────────────────
print(f"Status : {solver.getModelReturnStatus()}")
print(f"Optimal: {solver.getPrimalBound():.6f}  (known optimum ≈ 4.5796)")
print(f"\nPrimal solution improvements ({len(primal_log)} updates):")
print(f"  {'iter':>5}  {'objective':>12}  {'dual bound':>12}  {'gap %':>8}")
for row in primal_log:
    print(f"  {row['iter']:>5}  {row['obj']:>12.6f}  {row['dual']:>12.6f}  {row['gap (%)']:>8.4f}")

# ── Example 3: inspect an ExternalHyperplane object ──────────────────────
hp = SHOTpy.ExternalHyperplane()
hp.variableIndexes      = [0, 2]        # x1 and x3
hp.variableCoefficients = [1.5, -0.5]   # 1.5*x1 - 0.5*x3 <= rhs
hp.rhsValue             = 3.0
hp.isGlobal             = True
hp.source               = SHOTpy.HyperplaneSource.External
hp.description          = "my_custom_cut"

print("\nExample ExternalHyperplane:")
print(f"  indexes : {hp.variableIndexes}")
print(f"  coeffs  : {hp.variableCoefficients}")
print(f"  rhs     : {hp.rhsValue}")
print(f"  global  : {hp.isGlobal}")
print(f"  source  : {hp.source}")

# The callbacks were called, and SHOT found the known optimum. Since early_stop can terminate SHOT before the
# gap is closed, the optimum need not be proven, but the gap is below 1 %
assert solver.getModelReturnStatus() in (SHOTpy.ModelReturnStatus.OptimalGlobal,
                                          SHOTpy.ModelReturnStatus.FeasibleSolution)
assert solver.getRelativeObjectiveGap() < 0.01
assert abs(solver.getPrimalBound() - 4.579582) < 1e-3
assert len(primal_log) > 0

PrimalCandidateCheck callback

PrimalCandidateCheck is called for every primal candidate before SHOT checks whether it is feasible and better than the incumbent — unlike NewPrimalSolution, which is called only after a solution has been stored. Both contexts have these properties:

Property

Type

Description

iterationNumber

int

Current iteration

objectiveValue

float

Objective value of the candidate point

dualBound

float

Dual bound at this point in the search

relativeGap

float

Relative optimality gap

isMinimization

bool

True for minimization problems

point

list[float]

Candidate variable values

source

PrimalSolutionSource

Where the candidate comes from

Action: ctx.rejectCandidate() rejects the candidate outright (the feasibility check is skipped and the candidate is never recorded as a primal solution). A candidate that is not rejected is checked normally.

This callback is useful for:

  • Logging all candidate points examined by the primal heuristics

  • Filtering candidates based on problem-specific criteria — e.g. enforcing extra constraints that are not in the model (call ctx.rejectCandidate())

Rejecting candidates only affects the primal side: the dual problems SHOT solves still do not contain the extra constraints, so the dual bound is that of the model without them, and SHOT does not search for solutions that fulfil them. Rejection alone is therefore not a way to solve a model with constraints left out; the constraints that should hold should be in the model.

import SHOTpy

# ── Build the ex1223b MINLP (same problem used in the main callback example) ──
def _build():
    s = SHOTpy.Solver()
    s.updateSetting("Output.Console.LogLevel", 6)  # silent
    p = SHOTpy.Problem(s); p.name = "ex1223b_cand"
    x1 = p.addVariable("x1", SHOTpy.VariableType.Real, 0.0, 10.0)
    x2 = p.addVariable("x2", SHOTpy.VariableType.Real, 0.0, 10.0)
    x3 = p.addVariable("x3", SHOTpy.VariableType.Real, 0.0, 10.0)
    b4 = p.addVariable("b4", SHOTpy.VariableType.Binary)
    b5 = p.addVariable("b5", SHOTpy.VariableType.Binary)
    b6 = p.addVariable("b6", SHOTpy.VariableType.Binary)
    b7 = p.addVariable("b7", SHOTpy.VariableType.Binary)

    p.setObjective((b4 - 1.0)**2 + (b5 - 2.0)**2 + (b6 - 1.0)**2 - SHOTpy.log(1.0 + b7)
                     + (x1 - 1.0)**2 + (x2 - 2.0)**2 + (x3 - 3.0)**2, SHOTpy.ObjectiveDirection.Minimize)

    p.addConstraint(x1 + x2 + x3 + b4 + b5 + b6 <= 5, "e1")
    p.addConstraint(b6**2 + x1**2 + x2**2 + x3**2 <= 5.5, "e2")
    p.addConstraint(x1 + b4 <= 1.2, "e3")
    p.addConstraint(x2 + b5 <= 1.8, "e4")
    p.addConstraint(x3 + b6 <= 2.5, "e5")
    p.addConstraint(x1 + b7 <= 1.2, "e6")
    p.addConstraint(b5**2 + x2**2 <= 1.64, "e7")
    p.addConstraint(b6**2 + x3**2 <= 4.25, "e8")
    p.addConstraint(b5**2 + x3**2 <= 4.64, "e9")
    p.finalize(); s.setProblem(p)
    return s

# ── Run 1: accept all candidates (none is rejected) ───────────────────────────
s1 = _build()
cand1, acc1 = [], []
s1.registerCallback(SHOTpy.CallbackLocation.PrimalCandidateCheck,
                    lambda ctx: cand1.append(round(ctx.objectiveValue, 4)))
s1.registerCallback(SHOTpy.CallbackLocation.NewPrimalSolution,
                    lambda ctx: acc1.append(round(ctx.objectiveValue, 4)))
assert s1.solveProblem()
# ── Run 2: reject sub-optimal candidates (obj > 5) ────────────────────────────
s2 = _build()
cand2, acc2 = [], []

def on_candidate(ctx):
    cand2.append(round(ctx.objectiveValue, 4))
    # Skip the feasibility check for this candidate
    if ctx.objectiveValue > 5.0:
        ctx.rejectCandidate()

s2.registerCallback(SHOTpy.CallbackLocation.PrimalCandidateCheck, on_candidate)
s2.registerCallback(SHOTpy.CallbackLocation.NewPrimalSolution,
                    lambda ctx: acc2.append(round(ctx.objectiveValue, 4)))
assert s2.solveProblem()
# ── Compare results ────────────────────────────────────────────────────────────
print("── Accept all ────────────────────────────────────────────────")
print(f"  Candidates seen   : {cand1}")
print(f"  Incumbents found  : {acc1}")
print()
print("── Reject obj > 5 ────────────────────────────────────────────")
print(f"  Candidates seen   : {cand2}")
print(f"  Incumbents found  : {acc2}")
print(f"  (rejected {sum(1 for v in cand2 if v > 5.0)} candidate(s) with obj > 5)")

# Both runs find the optimum, but the second never accepts an incumbent with objective above 5
assert abs(s1.getPrimalBound() - 4.579582) < 1e-3
assert abs(s2.getPrimalBound() - 4.579582) < 1e-3
assert any(v > 5.0 for v in cand2)
assert all(v <= 5.0 for v in acc2)

Hyperplane Selection Callback: Solving shot_ex_jogo

This example demonstrates the HyperplaneSelection callback by solving a small MINLP whose problem data comes from test/data/shot_ex_jogo.gms. The problem is constructed directly through the Python API.

\[ \min_{x_1,\,x_2} \; -x_1 - x_2 \]

subject to

\[ 2x_1 - 3x_2 \le 2, \qquad x_1 \in [1,20],\; x_2 \in \mathbb{Z} \cap [1,20] \]
\[ 0.15(x_1-8)^2 + 0.1(x_2-6)^2 + \frac{0.025\,e^{x_1}}{x_2^2} \le 5 \]
\[ \frac{1}{x_1} + \frac{1}{x_2} - \sqrt{x_1}\,\sqrt{x_2} \le -4 \]

All built-in SHOT hyperplane generation is disabled (CutStrategy = OnlyExternal). The callback computes a linearisation cut (supporting hyperplane) for the most-violated nonlinear constraint at every iteration and adds it to SHOT as an ExternalHyperplane with ctx.addHyperplane().

import SHOTpy

# ── Create solver and configure settings ─────────────────────────────────────
solver = SHOTpy.Solver()

solver.updateSetting("Output.Console.LogLevel", 2)

# Disable all built-in hyperplane generation; cuts come exclusively from the callback
solver.updateSetting("Model.Convexity.AssumeConvex", True)   # True\n
solver.updateSetting("Model.Reformulation.Constraint.PartitionQuadraticTerms", 2)  # ES_PartitionNonlinearSums::Never
solver.updateSetting("Dual.CutStrategy", 2)    # ES_HyperplaneCutStrategy::OnlyExternal
solver.updateSetting("Dual.TreeStrategy", 0)    # ES_TreeStrategy::MultiTree

# ── Build the shot_ex_jogo problem via the Python API ─────────────────────────
#
#   min  -x1 - x2
#   s.t. 2*x1 - 3*x2                            <=  2   (l,  linear)
#        0.15*(x1-8)^2 + 0.1*(x2-6)^2
#          + 0.025*exp(x1)/x2^2                 <=  5   (c1, nonlinear)
#        1/x1 + 1/x2 - sqrt(x1)*sqrt(x2)        <= -4   (c2, nonlinear / signomial)
#   x1 in [1, 20]  (continuous)
#   x2 in [1, 20]  (integer)

problem      = SHOTpy.Problem(solver)
problem.name = "shot_ex_jogo"

x1 = problem.addVariable("x1", SHOTpy.VariableType.Real,    1.0, 20.0)
x2 = problem.addVariable("x2", SHOTpy.VariableType.Integer, 1.0, 20.0)

# Objective: minimize -x1 - x2
problem.setObjective(-x1 - x2, SHOTpy.ObjectiveDirection.Minimize)

# Constraint l: 2*x1 - 3*x2 <= 2
problem.addConstraint(2 * x1 - 3 * x2 <= 2, "l")

# Constraint c1: 0.15*(x1-8)^2 + 0.1*(x2-6)^2 + 0.025*exp(x1)/x2^2 <= 5
problem.addConstraint(0.15 * (x1 - 8.0)**2 + 0.1 * (x2 - 6.0)**2 + 0.025 * SHOTpy.exp(x1) / x2**2 <= 5, "c1")

# Constraint c2: 1/x1 + 1/x2 - x1^0.5 * x2^0.5 <= -4, which finalize() turns into signomial terms
problem.addConstraint(1 / x1 + 1 / x2 - x1**0.5 * x2**0.5 <= -4, "c2")

problem.finalize()

# Pass the same problem as both original and reformulated to bypass internal reformulation
assert solver.setProblem(problem, problem)
print("Problem created:")
print(problem.toString())

# ── Register external hyperplane callback ─────────────────────────────────────
#
# At every iteration SHOT calls this function with the solution points of the
# dual problem. For a point we linearise the most-violated nonlinear constraint
# at that point and add the resulting cut as an ExternalHyperplane.
#
# Supporting hyperplane for  g(x) <= b  at  x0:
#   grad g(x0) . x  <=  grad g(x0) . x0  -  (g(x0) - b)
#                   <=  sum_i grad_i * x0_i  -  violation

cut_log = []   # record (iteration, constraint name, violation) for display
iteration_log = {}  # callback points and cuts, grouped by SHOT iteration for plotting

def generate_hyperplanes(ctx):
    if not ctx.solutionPoints or ctx.iterationNumber == 0:
        return

    reform_problem = ctx.reformulatedProblem
    # Selection follows the preceding relaxation solve; its points belong to that solve iteration.
    solve_iteration = ctx.iterationNumber - 1
    snapshot = iteration_log.setdefault(solve_iteration, {"points": [], "cuts": []})

    for sol_point in ctx.solutionPoints:
        snapshot["points"].append(tuple(sol_point.point[:2]))
        dev_idx = sol_point.maxDeviation.index

        # Index must be valid inside the nonlinear-constraint list
        if dev_idx < 0 or dev_idx >= len(reform_problem.nonlinearConstraints):
            continue

        violation = sol_point.maxDeviation.value
        if violation <= 0.0:
            continue   # constraint is satisfied at this point

        # Looks up by the constraint's own .index attribute
        constraint = next(
            (nlc for nlc in reform_problem.nonlinearConstraints if nlc.index == dev_idx), None)
        
        point      = sol_point.point

        # Sparse gradient dict {var_index: coefficient}
        gradient = constraint.calculateGradient(point)

        if not gradient:
            continue

        var_indices  = list(gradient.keys())
        var_coeffs   = list(gradient.values())

        # rhs = grad . x0 - violation
        rhs = sum(gradient[i] * point[i] for i in var_indices) - violation

        hp                      = SHOTpy.ExternalHyperplane()
        hp.variableIndexes      = var_indices
        hp.variableCoefficients = var_coeffs
        hp.rhsValue             = rhs
        hp.isGlobal             = True
        hp.source               = SHOTpy.HyperplaneSource.External
        hp.description          = f"ext_hp_iter{ctx.iterationNumber}_c{dev_idx}"

        ctx.addHyperplane(hp)
        cut_log.append((ctx.iterationNumber, constraint.name, round(violation, 6)))
        snapshot["cuts"].append({
            "constraint": constraint.name,
            "point": tuple(point[:2]),
            "gradient": (gradient.get(0, 0.0), gradient.get(1, 0.0)),
            "rhs": rhs,
        })

        break   # one cut per callback invocation is sufficient

solver.registerCallback(SHOTpy.CallbackLocation.HyperplaneSelection, generate_hyperplanes)

# ── Solve ──────────────────────────────────────────────────────────────────────
print("\nSolving with external hyperplanes only...\n")
assert solver.solveProblem()
# ── Results ───────────────────────────────────────────────────────────────────
print(f"\nStatus : {solver.getModelReturnStatus()}")

if solver.getPrimalSolutions():
    sol = solver.getPrimalSolution()
    print(f"Objective value : {sol.objValue:.6f}")
    print(f"x1 = {sol.point[0]:.6f}  (continuous)")
    print(f"x2 = {sol.point[1]:.6f}  (integer)")
else:
    print("No primal solution found.")

print(f"\nExternal cuts added: {len(cut_log)}")
if cut_log:
    print(f"  {'iter':>6}  {'constraint':>12}  {'violation':>12}")
    for (it, cname, viol) in cut_log[:10]:
        print(f"  {it:>6}  {cname:>12}  {viol:>12.6f}")
    if len(cut_log) > 10:
        print(f"  ... ({len(cut_log) - 10} more)")

# The cuts of the callback alone give the known optimum
assert solver.getModelReturnStatus() == SHOTpy.ModelReturnStatus.OptimalGlobal
assert abs(solver.getPrimalSolution().objValue + 20.903615) < 1e-3
assert len(cut_log) > 0

The callback also records the relaxation solution points and the cuts it selects. The cell below draws one image for each SHOT iteration. The green region satisfies all three constraints with the integrality requirement on x2 relaxed; horizontal dotted lines mark integer x2 values. Thick green horizontal segments show the intersection of the relaxed feasible region with those integer rows: every point on a segment satisfies all three constraints and has integer x2. Dashed blue curves are contours of the objective -x1-x2. Black, purple, and orange curves are the original constraint boundaries. Gray points and faint red lines come from earlier iterations; the current callback points and newly selected cuts are emphasized. Orange stars show the linearization points. The selection callback runs after a relaxation solve, so its points and newly selected cuts are shown in that solve iteration’s image; SHOT applies those cuts to subsequent relaxations. The final image also marks the accepted primal solution. The axes remain fixed across images for comparison.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from matplotlib.patches import Patch

# Integer-relaxed feasible region on the two original variables.
x_grid = np.linspace(1.0, 20.0, 240)
y_grid = np.linspace(1.0, 20.0, 240)
X, Y = np.meshgrid(x_grid, y_grid)
linear = 2 * X - 3 * Y - 2
c1 = 0.15 * (X - 8)**2 + 0.1 * (Y - 6)**2 + 0.025 * np.exp(X) / Y**2 - 5
c2 = 1 / X + 1 / Y - np.sqrt(X * Y) + 4
feasible = (linear <= 0) & (c1 <= 0) & (c2 <= 0)
objective = -X - Y

# Feasible portions of each integer-x2 row. x1 remains continuous, so these are segments.
integer_x = np.linspace(1.0, 20.0, 1500)
integer_feasible_rows = []
for integer_y in range(1, 21):
    row_linear = 2 * integer_x - 3 * integer_y <= 2
    row_c1 = (0.15 * (integer_x - 8)**2 + 0.1 * (integer_y - 6)**2
              + 0.025 * np.exp(integer_x) / integer_y**2 <= 5)
    row_c2 = (1 / integer_x + 1 / integer_y - np.sqrt(integer_x * integer_y) <= -4)
    row_x = np.where(row_linear & row_c1 & row_c2, integer_x, np.nan)
    if np.isfinite(row_x).any():
        integer_feasible_rows.append((integer_y, row_x))

number_of_iterations = solver.getSolutionStatistics().numberOfIterations
for iteration in range(1, number_of_iterations + 1):
    fig, ax = plt.subplots(figsize=(8.5, 5.1), constrained_layout=True)
    ax.contourf(X, Y, feasible.astype(int), levels=[-0.5, 0.5, 1.5],
                colors=["#f8f1ee", "#d7efe1"])
    objective_lines = ax.contour(X, Y, objective, levels=np.arange(-36, -3, 4),
                                 colors="#6294c5", linestyles="dashed", linewidths=0.8)
    ax.clabel(objective_lines, inline=True, fontsize=7, fmt="%g")
    ax.contour(X, Y, linear, levels=[0], colors="#333333", linewidths=1.2)
    ax.contour(X, Y, c1, levels=[0], colors="#d78a27", linewidths=1.2)
    ax.contour(X, Y, c2, levels=[0], colors="#8856a7", linewidths=1.2)
    for integer_value in range(1, 21):
        ax.axhline(integer_value, color="#777777", linestyle=":", linewidth=0.35, alpha=0.35)
    for integer_y, row_x in integer_feasible_rows:
        ax.plot(row_x, np.full_like(row_x, integer_y), color="#16834e", linewidth=3,
                solid_capstyle="round", zorder=3)

    for recorded_iteration in sorted(iteration_log):
        if recorded_iteration > iteration:
            break
        snapshot = iteration_log[recorded_iteration]
        current = recorded_iteration == iteration
        for cut in snapshot["cuts"]:
            gx, gy = cut["gradient"]
            ax.contour(X, Y, gx * X + gy * Y - cut["rhs"], levels=[0],
                       colors=["#d62728" if current else "#b98c8c"],
                       linewidths=1.8 if current else 0.7, alpha=0.95 if current else 0.4)
            px, py = cut["point"]
            ax.scatter(px, py, marker="*", s=125 if current else 35,
                       color="#ee7b19" if current else "#c7a58a", zorder=5)
        for px, py in snapshot["points"]:
            ax.scatter(px, py, facecolors="none", edgecolors="black" if current else "#777777",
                       s=75 if current else 25, linewidths=1.4 if current else 0.7, zorder=6)

    if iteration == number_of_iterations and solver.getPrimalSolutions():
        final_point = solver.getPrimalSolution().point
        ax.scatter(final_point[0], final_point[1], marker="D", s=65,
                   color="#159447", edgecolors="black", linewidths=0.5, zorder=7)

    ax.set(xlim=(0.8, 20.4), ylim=(0.8, 20.4), xlabel="x1 (continuous)", ylabel="x2 (integer in the MINLP)",
           title=f"SHOT iteration {iteration}/{number_of_iterations}: "
                 f"{sum(len(v['cuts']) for k, v in iteration_log.items() if k <= iteration)} cuts selected")
    ax.set_aspect("equal", adjustable="box")
    ax.legend(handles=[
        Patch(facecolor="#d7efe1", label="Relaxed feasible region"),
        Line2D([], [], color="#16834e", linewidth=3, label="Integer-feasible segments"),
        Line2D([], [], color="#6294c5", linestyle="--", label="Objective contours"),
        Line2D([], [], color="#333333", label="Linear constraint"),
        Line2D([], [], color="#d78a27", label="c1 boundary"),
        Line2D([], [], color="#8856a7", label="c2 boundary"),
        Line2D([], [], marker="o", markerfacecolor="none", color="black", linestyle="None",
               label="Current relaxation point"),
        Line2D([], [], marker="*", color="#ee7b19", linestyle="None", label="Linearization point"),
        Line2D([], [], color="#d62728", label="New cut"),
        Line2D([], [], color="#b98c8c", label="Earlier cuts"),
        Line2D([], [], marker="D", color="#159447", linestyle="None", label="Final primal solution"),
    ], loc="upper left", bbox_to_anchor=(1.02, 1), fontsize=7)
    plt.show()
    plt.close(fig)

Dual Bound and Primal Solution Callbacks

Two more callbacks let the user feed external information directly into the search:

Location

Action

Effect

DualBoundUpdate

ctx.setDualBound(value)

Propose a dual bound; SHOT uses it if it is better than its own.

PrimalCandidateSearch

ctx.addPrimalSolution(point)

Add a candidate primal point; can be called several times.

Both contexts have the common properties, e.g., iterationNumber, dualBound and primalBound.

Linear mixed-integer problems honor Dual.TreeStrategy with every backend (Cbc and HiGHS support multi-tree). The direct MIQCQP strategy is used only for quadratic reformulated problems; both callbacks run before and after the CPLEX or Gurobi solve in that strategy. A feasible point added before the solve becomes a primal solution; SHOT can then pass it as a MIP start and use its objective as a cutoff. A point added after the solve cannot affect that completed MIP search. An external dual bound updates SHOT’s gap and termination checks; it does not change CPLEX’s or Gurobi’s internal search bound.

The example below solves shot_ex_jogo again, this time with both callbacks active:

  1. DualBoundUpdate – waits until iteration 3 before proposing the known optimal dual bound −20.9036, tightening the gap early.

  2. PrimalCandidateSearch – on the very first call adds the known optimal point [8.903615, 12.0] (previously verified by a warm-up solve), letting the solver skip most of the primal search.

Together they let the solver terminate in iteration 3, when the dual bound is proposed, since the known optimal solution has already been added in iteration 1. A proposed dual bound must be valid, i.e., no feasible solution may be better than it, since SHOT uses it to decide that the solution is optimal; the same holds for cuts added in the callbacks, which must not cut off any feasible solution.

import SHOTpy
import math

# ── Warm-up solve to collect a verified feasible primal solution ──────────────
#
# We first solve the problem normally to get a solution point that SHOT has
# already checked for feasibility.  This avoids floating-point rounding issues
# when re-injecting the point in the main solve.

ws = SHOTpy.Solver()
ws.updateSetting("Output.Console.LogLevel", 6)   # silent
ws.updateSetting("Model.Convexity.AssumeConvex", True)

def build_shot_ex_jogo(solver):
    problem = SHOTpy.Problem(solver)
    problem.name = "shot_ex_jogo"

    x1 = problem.addVariable("x1", SHOTpy.VariableType.Real,    1.0, 20.0)
    x2 = problem.addVariable("x2", SHOTpy.VariableType.Integer, 1.0, 20.0)

    problem.setObjective(-x1 - x2, SHOTpy.ObjectiveDirection.Minimize)

    problem.addConstraint(2 * x1 - 3 * x2 <= 2, "l")

    problem.addConstraint(0.15 * (x1 - 8.0)**2 + 0.1 * (x2 - 6.0)**2 + 0.025 * SHOTpy.exp(x1) / x2**2 <= 5, "c1")

    problem.addConstraint(1 / x1 + 1 / x2 - x1**0.5 * x2**0.5 <= -4, "c2")

    problem.finalize()
    return problem

assert ws.setProblem(build_shot_ex_jogo(ws), build_shot_ex_jogo(ws))
assert ws.solveProblem()
known_point = list(ws.getPrimalSolution().point)   # verified-feasible optimal point
known_obj   = ws.getPrimalSolution().objValue
known_dual  = -20.903647306104258                  # from reference solution

print(f"Warm-up solve finished.")
print(f"  Optimal point : x1={known_point[0]:.6f}, x2={known_point[1]:.1f}")
print(f"  Optimal obj   : {known_obj:.6f}")
print(f"  Known dual    : {known_dual:.6f}")

# ── Main solve with DualBoundUpdate and PrimalCandidateSearch callbacks ──────

solver = SHOTpy.Solver()
solver.updateSetting("Output.Console.LogLevel", 2)
solver.updateSetting("Model.Convexity.AssumeConvex", True)
solver.updateSetting("Dual.Relaxation.Use", False)
solver.updateSetting("Dual.TreeStrategy", 0)    # MultiTree
assert solver.setProblem(build_shot_ex_jogo(solver), build_shot_ex_jogo(solver))
# ── Callback state ────────────────────────────────────────────────────────────
INJECT_DUAL_AFTER_ITER = 3   # wait this many iterations before tightening the dual bound

dual_log    = []
primal_log  = []
event_log   = []

# ── DualBoundUpdate callback ──────────────────────────────────────────────────
#
# Called once per main iteration.  We wait until iteration > INJECT_DUAL_AFTER_ITER
# so that the solver has already made some progress on its own, then propose the
# known optimal dual bound to close the gap immediately.

dual_injected = [False]

def provide_dual_bound(ctx):
    entry = {
        "iter":    ctx.iterationNumber,
        "current": ctx.dualBound,
        "injected": None,
    }
    if (not dual_injected[0]
            and ctx.iterationNumber >= INJECT_DUAL_AFTER_ITER
            and ctx.dualBound < known_dual):
        dual_injected[0]  = True
        entry["injected"] = known_dual
        event_log.append(f"iter {ctx.iterationNumber:>2}: [DualBoundUpdate]  "
                         f"current={ctx.dualBound:.4f}  "
                         f"→ proposing known dual {known_dual:.6f}")
        dual_log.append(entry)
        ctx.setDualBound(known_dual)
        return

    event_log.append(f"iter {ctx.iterationNumber:>2}: [DualBoundUpdate]  "
                     f"current={ctx.dualBound:.4f}  → (no update)")
    dual_log.append(entry)   # no bound is proposed

# ── PrimalCandidateSearch callback ────────────────────────────────────────────
#
# Called once per main iteration.  On the very first call we add the
# warm-up optimal point; SHOT validates it and, if feasible, accepts it as the
# new incumbent (and calls the NewPrimalSolution callback).

primal_injected = [False]

def provide_primal_solution(ctx):
    entry = {
        "iter":    ctx.iterationNumber,
        "current": ctx.primalBound,
        "injected": None,
    }
    if not primal_injected[0]:
        primal_injected[0]  = True
        entry["injected"]   = known_point
        event_log.append(f"iter {ctx.iterationNumber:>2}: [PrimalCandidateSearch]  "
                         f"current={ctx.primalBound:.4f}  "
                         f"→ adding known point {[round(v,4) for v in known_point]}")
        primal_log.append(entry)
        ctx.addPrimalSolution(known_point)   # list[float] with one value per variable
        return

    event_log.append(f"iter {ctx.iterationNumber:>2}: [PrimalCandidateSearch]  "
                     f"current={ctx.primalBound:.4f}  → (no update)")
    primal_log.append(entry)

# ── NewPrimalSolution notification ───────────────────────────────────────────
def on_new_primal(ctx):
    event_log.append(f"iter {ctx.iterationNumber:>2}: [NewPrimalSolution]  "
                     f"obj={ctx.objectiveValue:.6f}  "
                     f"point={[round(v,4) for v in ctx.point]}")

L = SHOTpy.CallbackLocation
solver.registerCallback(L.NewPrimalSolution,     on_new_primal)
solver.registerCallback(L.DualBoundUpdate,       provide_dual_bound)
solver.registerCallback(L.PrimalCandidateSearch, provide_primal_solution)

print(f"\nSolving with DualBoundUpdate (propose after iter {INJECT_DUAL_AFTER_ITER}) "
      f"and PrimalCandidateSearch (add on first call)...\n")
assert solver.solveProblem()
# ── Print event log ───────────────────────────────────────────────────────────
print("\nCallback events (in order):")
for msg in event_log:
    print(f"  {msg}")

# ── Final results ─────────────────────────────────────────────────────────────
print(f"\nStatus    : {solver.getModelReturnStatus()}")
if solver.getPrimalSolutions():
    sol = solver.getPrimalSolution()
    print(f"Objective : {sol.objValue:.6f}  (known optimum {known_obj:.6f})")
    print(f"x1 = {sol.point[0]:.6f}  x2 = {sol.point[1]:.1f}")
    print(f"Iterations: {solver.getSolutionStatistics().numberOfIterations}")
else:
    print("No primal solution found.")

# The dual bound was proposed, and the solution is the known optimum
assert dual_injected[0]
assert solver.getModelReturnStatus() == SHOTpy.ModelReturnStatus.OptimalGlobal
assert abs(solver.getPrimalSolution().objValue - known_obj) < 1e-4

Plotting primal and dual bounds over time

BoundsUpdated reports changes to SHOT’s recorded bounds. The callback below collects every update and refreshes one notebook figure about twice a second while SHOT is solving. The final refresh includes all events, even those between screen updates. The title shows the current absolute gap and relative gap (as a percentage), calculated from the current dual and primal bounds. The plot uses the fo7.osil benchmark bundled with this repository, so run the notebook from the repository or one of its build directories. Install matplotlib if needed.

This example uses one MIP thread to make the search easier to watch and limits the solve to 35 seconds. On this machine, the debug build took about 20 seconds and produced thousands of bound updates; run time varies with hardware, solver, settings, and build type. A 35-second time limit is a maximum, not a guaranteed duration. The time axis expands automatically if a run continues past 35 seconds.

import math
import time
from pathlib import Path
from IPython.display import display
import matplotlib.pyplot as plt
import SHOTpy

# Find the bundled benchmark from the repository root or a build directory.
relative_file = Path("test/data/instances/MINLP-convex/fo7.osil")
problem_file = next((base / relative_file for base in (Path.cwd(), *Path.cwd().parents)
                     if (base / relative_file).is_file()), None)
if problem_file is None:
    raise FileNotFoundError(f"Run this notebook from the SHOT repository or a build directory: {relative_file}")

solver = SHOTpy.Solver()
solver.updateSetting("Termination.TimeLimit", 35.0)
solver.updateSetting("Dual.MIP.NumberOfThreads", 1)
solver.updateSetting("Output.Console.LogLevel", 6)
if not solver.setProblem(str(problem_file)):
    raise RuntimeError(f"Could not load {problem_file}")

samples = []
fig, ax = plt.subplots(figsize=(9, 5))
ax.set_xlabel("Elapsed solve time (s)")
ax.set_ylabel("Objective value")
ax.set_xlim(0, 35)
plot_handle = display(fig, display_id=True)
plt.close(fig)  # The display handle owns the visible notebook output.
last_refresh = time.monotonic()

def refresh_plot():
    ax.clear()
    if samples:
        times, dual, primal, absolute_gaps, relative_gaps = zip(*samples)
        ax.step(times, [v if math.isfinite(v) else math.nan for v in dual],
                where="post", label="Dual bound")
        ax.step(times, [v if math.isfinite(v) else math.nan for v in primal],
                where="post", label="Primal bound")
        ax.legend()
    ax.set_xlim(0, max(35.0, samples[-1][0] * 1.05) if samples else 35.0)
    ax.set_xlabel("Elapsed solve time (s)")
    ax.set_ylabel("Objective value")
    if samples:
        absolute_gap = f"{absolute_gaps[-1]:.6g}" if math.isfinite(absolute_gaps[-1]) else "unavailable"
        relative_gap = f"{100 * relative_gaps[-1]:.3g}%" if math.isfinite(relative_gaps[-1]) else "unavailable"
    else:
        absolute_gap = relative_gap = "unavailable"
    ax.set_title(f"SHOT bounds ({len(samples) - 1 if samples else 0} updates) | "
                 f"absolute gap: {absolute_gap} | relative gap: {relative_gap}")
    if plot_handle is not None:  # Outside Jupyter, display() has no update handle.
        plot_handle.update(fig)

def record_bounds(ctx: SHOTpy.BoundsUpdatedContext):
    global last_refresh
    if not samples:
        samples.append((0.0, ctx.previousDualBound, ctx.previousPrimalBound, math.inf, math.inf))
    samples.append((ctx.elapsedTime, ctx.dualBound, ctx.primalBound, ctx.absoluteGap, ctx.relativeGap))
    now = time.monotonic()
    if now - last_refresh >= 0.5:
        refresh_plot()
        last_refresh = now

solver.registerCallback(SHOTpy.CallbackLocation.BoundsUpdated, record_bounds)
try:
    solver.solveProblem()
finally:
    refresh_plot()  # Show the last updates even if the solve ends between refreshes.
print(f"Recorded {max(0, len(samples) - 1)} bound updates.")

ESH Interior Point Callback: Injecting Custom Interior Points

The InteriorPointSearch callback location lets you control which interior points are used by the Extended Supporting Hyperplane (ESH) algorithm. The ESH method needs a strictly interior point of the feasible region to generate supporting hyperplanes; normally SHOT finds this point automatically by solving a cutting-plane minimax NLP. The callback gives you two complementary powers:

  1. Inspect / filter the point(s) found by the internal strategy — the callback is called with ctx.interiorPoints containing the points the solver just found. Leave them unchanged by doing nothing, or replace them with ctx.setInteriorPoints(points). If the search found no usable point, this list is empty and you can supply one.

  2. Bypass the internal search entirely by setting Dual.ESH.InteriorPoint.Strategy = 1 (“Only external (through callback)”). The solver will skip the minimax NLP and rely completely on the points you set in the callback. This is useful when you already know a feasible interior point — for example from a previous solve.

Dual.ESH.InteriorPoint.Strategy
  0  Use internal strategy (default) — callback fires after the search, even if it finds no point
  1  Only external (through callback) — internal NLP is skipped entirely

Callback signature:

def on_interior_points(ctx: SHOTpy.InteriorPointSearchContext) -> None:
    ...

Property / action

Type

Description

ctx.interiorPoints

list[list[float]]

Current interior points in reformulated variables; empty if none were found or OnlyExternal is used

ctx.setInteriorPoints(points)

Replaces the points; each point can be in the variables of the original or the reformulated problem

ctx.originalProblem

Problem

The original problem object

ctx.reformulatedProblem

Problem

The reformulated problem used internally by SHOT

ctx.solutionStatistics

SolutionStatistics

Current solver statistics

The example below uses a two-phase pattern:

  • Phase 1 — solve ex1223b normally, capture the interior point via the callback.

  • Phase 2 — create a fresh solver with Dual.ESH.InteriorPoint.Strategy = OnlyExternal, inject the captured point, and verify the same optimal solution is reached.

import SHOTpy

# ── Helper: build ex1223b ─────────────────────────────────────────────────
def build_ex1223b_nb(solver):
    problem = SHOTpy.Problem(solver)
    problem.name = "ex1223b"
    x1 = problem.addVariable("x1", SHOTpy.VariableType.Real, 0.0, 10.0)
    x2 = problem.addVariable("x2", SHOTpy.VariableType.Real, 0.0, 10.0)
    x3 = problem.addVariable("x3", SHOTpy.VariableType.Real, 0.0, 10.0)
    b4 = problem.addVariable("b4", SHOTpy.VariableType.Binary)
    b5 = problem.addVariable("b5", SHOTpy.VariableType.Binary)
    b6 = problem.addVariable("b6", SHOTpy.VariableType.Binary)
    b7 = problem.addVariable("b7", SHOTpy.VariableType.Binary)

    problem.setObjective((b4 - 1.0)**2 + (b5 - 2.0)**2 + (b6 - 1.0)**2 - SHOTpy.log(1.0 + b7)
                     + (x1 - 1.0)**2 + (x2 - 2.0)**2 + (x3 - 3.0)**2, SHOTpy.ObjectiveDirection.Minimize)

    problem.addConstraint(x1 + x2 + x3 + b4 + b5 + b6 <= 5, "e1")
    problem.addConstraint(b6**2 + x1**2 + x2**2 + x3**2 <= 5.5, "e2")
    problem.addConstraint(x1 + b4 <= 1.2, "e3")
    problem.addConstraint(x2 + b5 <= 1.8, "e4")
    problem.addConstraint(x3 + b6 <= 2.5, "e5")
    problem.addConstraint(x1 + b7 <= 1.2, "e6")
    problem.addConstraint(b5**2 + x2**2 <= 1.64, "e7")
    problem.addConstraint(b6**2 + x3**2 <= 4.25, "e8")
    problem.addConstraint(b5**2 + x3**2 <= 4.64, "e9")
    problem.finalize()
    return problem

# The constraints of ex1223b are quadratic, which the MIP solver would otherwise handle directly, leaving ESH no
# nonlinear constraints and no interior point to search for. They are therefore treated as nonlinear.

# ═══════════════════════════════════════════════════════════════════════════
# Phase 1 — Solve normally, capture the interior point via callback
# ═══════════════════════════════════════════════════════════════════════════
print("Phase 1: solving ex1223b with default ESH strategy")
solver1 = SHOTpy.Solver()
solver1.updateSetting("Output.Console.LogLevel", 2)
solver1.updateSetting("Model.Convexity.AssumeConvex", True)
solver1.updateSetting("Model.Reformulation.Quadratics.Strategy", 0)  # Nonlinear, see below
assert solver1.setProblem(build_ex1223b_nb(solver1))
captured_points = []

def capture_cb(ctx):
    """Called after internal strategy finds interior point(s); we just record them."""
    captured_points.extend(ctx.interiorPoints)
    print(f"  [Phase 1 callback] {len(ctx.interiorPoints)} interior point(s) found internally")
    # No points are set, so they are kept unchanged

solver1.registerCallback(SHOTpy.CallbackLocation.InteriorPointSearch, capture_cb)
assert solver1.solveProblem()
phase1_obj = solver1.getPrimalSolution().objValue
print(f"Phase 1 optimal objective: {phase1_obj:.6f}  (expected ≈ 4.5796)")
print(f"Captured {len(captured_points)} interior point(s) to use in Phase 2")

# ═══════════════════════════════════════════════════════════════════════════
# Phase 2 — Skip internal NLP; inject the captured point via callback
# ═══════════════════════════════════════════════════════════════════════════
print("\nPhase 2: solving with Dual.ESH.InteriorPoint.Strategy = OnlyExternal")
solver2 = SHOTpy.Solver()
solver2.updateSetting("Output.Console.LogLevel", 2)
solver2.updateSetting("Model.Convexity.AssumeConvex", True)
solver2.updateSetting("Model.Reformulation.Quadratics.Strategy", 0)  # Nonlinear, see below
solver2.updateSetting("Dual.ESH.InteriorPoint.Strategy", 1)  # OnlyExternal
assert solver2.setProblem(build_ex1223b_nb(solver2))
def inject_cb(ctx):
    """Provide the captured Phase-1 interior point; internal NLP was skipped."""
    print(f"  [Phase 2 callback] interiorPoints={len(ctx.interiorPoints)} "
          f"(should be 0 for OnlyExternal)")
    ctx.setInteriorPoints(captured_points)  # inject the captured points

solver2.registerCallback(SHOTpy.CallbackLocation.InteriorPointSearch, inject_cb)
assert solver2.solveProblem()
phase2_obj = solver2.getPrimalSolution().objValue
print(f"Phase 2 optimal objective: {phase2_obj:.6f}")
print(f"Objectives match: {abs(phase2_obj - phase1_obj) < 1e-4}")

# The interior points from Phase 1 give the same optimum in Phase 2
assert len(captured_points) > 0
assert abs(phase1_obj - 4.579582) < 1e-3
assert abs(phase2_obj - phase1_obj) < 1e-4

ESH Interior Point via Auxiliary Problem

When no interior point is available from a previous solve — for example at the very start of a workflow — you can find one by solving a small auxiliary NLP.

Idea: take the original problem, relax all integer/binary variables to continuous, add an auxiliary slack variable \(\mu\), subtract \(\mu\) from every nonlinear constraint, and minimise \(\mu\). If the optimal \(\mu^\star < 0\) the solution point is strictly interior to all nonlinear constraints, and can be fed directly into the InteriorPointSearch callback.

Modified problem (auxiliary)
──────────────────────────────────────────────────────────────────────
min  μ
s.t.  g_i(x) − μ  ≤  b_i   for each nonlinear constraint i
      all original linear constraints
      binary variables relaxed to Real [0, 1]
      μ ∈ [−100, 100]
──────────────────────────────────────────────────────────────────────

The cell below applies this to ex1223b:

  1. Auxiliary solve — builds the modified problem, solves it, and extracts the 7-variable interior point by dropping \(\mu\) (the last variable).

  2. Main solve — re-solves the original ex1223b with Dual.ESH.InteriorPoint.Strategy = OnlyExternal, injecting the auxiliary-problem interior point through the callback.

import SHOTpy

# ── Helper reused from above ─────────────────────────────────────────────────
# build_ex1223b_nb(solver) is defined in the earlier ESH cell; reuse it here.

# ══════════════════════════════════════════════════════════════════════════════
# Phase 1 — Auxiliary problem: minimise mu to find a strictly interior point
# ══════════════════════════════════════════════════════════════════════════════
print("Phase 1: solving auxiliary minimize-mu problem to find an interior point")

solver_aux = SHOTpy.Solver()
solver_aux.updateSetting("Output.Console.LogLevel", 6)  # silent
solver_aux.updateSetting("Model.Reformulation.Quadratics.Strategy", 0)  # Nonlinear

problem_aux = SHOTpy.Problem(solver_aux)
problem_aux.name = "ex1223b_interior"

# Original 7 variables — binary b4-b7 relaxed to Real [0, 1]
x1 = problem_aux.addVariable("x1", SHOTpy.VariableType.Real, 0.0, 10.0)
x2 = problem_aux.addVariable("x2", SHOTpy.VariableType.Real, 0.0, 10.0)
x3 = problem_aux.addVariable("x3", SHOTpy.VariableType.Real, 0.0, 10.0)
b4 = problem_aux.addVariable("b4", SHOTpy.VariableType.Real, 0.0, 1.0)
b5 = problem_aux.addVariable("b5", SHOTpy.VariableType.Real, 0.0, 1.0)
b6 = problem_aux.addVariable("b6", SHOTpy.VariableType.Real, 0.0, 1.0)
b7 = problem_aux.addVariable("b7", SHOTpy.VariableType.Real, 0.0, 1.0)
# Auxiliary slack variable
mu = problem_aux.addVariable("mu", SHOTpy.VariableType.Real, -100.0, 100.0)

# Objective: minimize mu
problem_aux.setObjective(mu, SHOTpy.ObjectiveDirection.Minimize)

# e1: linear, unchanged
problem_aux.addConstraint(x1 + x2 + x3 + b4 + b5 + b6 <= 5, "e1")

# e2, e7, e8, e9: quadratic constraints with -mu subtracted
problem_aux.addConstraint(b6**2 + x1**2 + x2**2 + x3**2 - mu <= 5.5, "e2")

# e3-e6: linear, unchanged
problem_aux.addConstraint(x1 + b4 <= 1.2, "e3")
problem_aux.addConstraint(x2 + b5 <= 1.8, "e4")
problem_aux.addConstraint(x3 + b6 <= 2.5, "e5")
problem_aux.addConstraint(x1 + b7 <= 1.2, "e6")

problem_aux.addConstraint(b5**2 + x2**2 - mu <= 1.64, "e7")
problem_aux.addConstraint(b6**2 + x3**2 - mu <= 4.25, "e8")
problem_aux.addConstraint(b5**2 + x3**2 - mu <= 4.64, "e9")

problem_aux.finalize()
assert solver_aux.setProblem(problem_aux)
assert solver_aux.solveProblem()
aux_sol      = solver_aux.getPrimalSolution()
optimal_mu   = aux_sol.objValue
interior_pt  = list(aux_sol.point[:7])   # drop mu (index 7)

print(f"Optimal mu : {optimal_mu:.6f}")
if optimal_mu < 0:
    print(f"  mu < 0  → strictly interior to all nonlinear constraints ✓")
else:
    print(f"  Warning: mu ≥ 0 — point may not be strictly interior")
print(f"Interior point (7 vars): {[round(v, 4) for v in interior_pt]}")

# ══════════════════════════════════════════════════════════════════════════════
# Phase 2 — Solve original ex1223b with OnlyExternal, injecting the aux point
# ══════════════════════════════════════════════════════════════════════════════
print("\nPhase 2: solving ex1223b with Dual.ESH.InteriorPoint.Strategy = OnlyExternal")
print("         injecting the auxiliary-problem interior point via callback")

solver2 = SHOTpy.Solver()
solver2.updateSetting("Output.Console.LogLevel", 6)
solver2.updateSetting("Model.Reformulation.Quadratics.Strategy", 0)  # Nonlinear
solver2.updateSetting("Dual.ESH.InteriorPoint.Strategy", 1)          # OnlyExternal
assert solver2.setProblem(build_ex1223b_nb(solver2))
def inject_aux_point(ctx):
    print(f"  [Callback] called — interiorPoints={len(ctx.interiorPoints)} "
          f"(should be 0 for OnlyExternal)")
    ctx.setInteriorPoints([interior_pt])   # one interior point

solver2.registerCallback(SHOTpy.CallbackLocation.InteriorPointSearch, inject_aux_point)
assert solver2.solveProblem()
phase2_obj = solver2.getPrimalSolution().objValue
print(f"\nPhase 2 objective  : {phase2_obj:.6f}  (known optimum ≈ 4.5796)")
print(f"Objectives match   : {abs(phase2_obj - 4.579582) < 0.01}")

# The point of the auxiliary problem is strictly interior, and gives the known optimum
assert optimal_mu < 0
assert abs(phase2_obj - 4.579582) < 1e-3