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 the interior point search of the ESH strategy

setInteriorPoints(points)

DualBoundUpdate

DualBoundUpdateContext

When SHOT can accept a dual bound

setDualBound(value)

PrimalCandidateSearch

PrimalCandidateSearchContext

When SHOT collects primal candidates

addPrimalSolution(point)

PrimalCandidateCheck

PrimalCandidateCheckContext

Before a primal candidate is checked

rejectCandidate()

NewPrimalSolution

NewPrimalSolutionContext

After a primal solution has been stored

–

TerminationCheck

TerminationCheckContext

When SHOT checks whether to terminate

–

HyperplaneSelection

HyperplaneSelectionContext

When SHOT selects the hyperplanes to add

addHyperplane(hyperplane)

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

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

    reform_problem = ctx.reformulatedProblem

    for sol_point in ctx.solutionPoints:
        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)))

        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

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.

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

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).

  2. Bypass the internal search entirely by setting 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.

ESH.InteriorPoint.Strategy
  0  Use internal strategy (default) — callback fires after, allowing modification
  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]]

Points found by the internal strategy (empty if OnlyExternal), in the variables of the reformulated problem

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 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 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 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 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