Continuous Optimization: Nonlinear Programming

# This code cell installs packages on Colab

import sys

if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper

    helper.easy_install()
else:
    sys.path.insert(0, "../")
    import helper
helper.set_plotting_style()

# `helper` is the ONE course module: the Colab install above, the house figure
# style, and the extract / archive / save-the-figure plumbing this notebook uses
# further down. It used to take a second `wget` for `pyomo_results.py`; that
# module is now a shim over this one, and one file is one thing that can fail.

import random

import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import pyomo.environ as pyo

# Seed the random number generator so this notebook is reproducible
# (Pyomo style guide, section 8).
random.seed(0)

def create_circle_model(circle_radii):
    """Create circle optimization model in Pyomo

    Arguments:
        circle_radii: dictionary with keys=circle name and value=radius (float)

    Returns:
        model: Pyomo model
    """
    model = pyo.ConcreteModel()
    # Ordered circle identifiers; radii and coordinates use metres.
    model.CIRCLES = pyo.Set(initialize=circle_radii.keys())
    model.R = pyo.Param(model.CIRCLES, domain=pyo.PositiveReals,
                        initialize=circle_radii, units=pyo.units.m)
    model.box_height = pyo.Var(domain=pyo.PositiveReals, units=pyo.units.m)
    model.box_width = pyo.Var(domain=pyo.PositiveReals, units=pyo.units.m)
    model.x = pyo.Var(model.CIRCLES, domain=pyo.PositiveReals, units=pyo.units.m)
    model.y = pyo.Var(model.CIRCLES, domain=pyo.PositiveReals, units=pyo.units.m)
    model.obj = pyo.Objective(expr=2 * (model.box_height + model.box_width),
                          sense=pyo.minimize)

    # One indexed rule creates one row per circle; b is the model.
    @model.Constraint(model.CIRCLES)
    def left_x_con(b, c):
        """Keep circle c inside the left edge of the box."""
        return b.x[c] >= b.R[c]

    @model.Constraint(model.CIRCLES)
    def left_y_con(b, c):
        """Keep circle c above the bottom edge of the box."""
        return b.y[c] >= b.R[c]

    @model.Constraint(model.CIRCLES)
    def right_x_con(b, c):
        """Keep circle c inside the right edge of the box."""
        return b.x[c] <= b.box_width - b.R[c]

    @model.Constraint(model.CIRCLES)
    def right_y_con(b, c):
        """Keep circle c below the top edge of the box."""
        return b.y[c] <= b.box_height - b.R[c]

    # Skip self-pairs and reversed duplicate pairs.
    @model.Constraint(model.CIRCLES, model.CIRCLES)
    def no_overlap_con(b, c1, c2):
        """Separate each distinct pair of circles."""
        if c1 < c2:
            return ((b.x[c1] - b.x[c2]) ** 2 + (b.y[c1] - b.y[c2]) ** 2
                    >= (b.R[c1] + b.R[c2]) ** 2)
        else:
            return pyo.Constraint.Skip
    return model

def initialize_circle_model(model, height_init=25, width_init=25, spread=10):
    """Initialize the x and y coordinates using uniform distribution

    Arguments:
        model: Pyomo model
        height_init: initial value for box_height (default=25)
        width_init: initial value for box_width (default=25)
        spread: width of the uniform window the centers are drawn from
            (default=10)

    Returns:
        Nothing. But per Pyomo scoping rules, the input argument `model`
        can be modified in this function.

    """
    # Initialize
    model.box_height = height_init
    model.box_width = width_init

    for i in model.CIRCLES:
        # Adding circle radii ensures the circle remains in the >0, >0 quadrant
        model.x[i] = random.uniform(0, spread) + pyo.value(model.R[i])
        model.y[i] = random.uniform(0, spread) + pyo.value(model.R[i])

# Create dictionary with circle data
circle_data = {"A": 10.0, "B": 5.0, "C": 3.0}
circle_data

# Access the keys
circle_data.keys()

# Create model
model = create_circle_model(circle_data)
model.pprint()

from pyomo.util.check_units import assert_units_consistent

# Raises InconsistentUnitsError if any constraint or the objective mixes units
assert_units_consistent(model)
print("Units are consistent.")

# Initialize model
initialize_circle_model(model)
model.pprint()

# ---------- EXTRACT, then PLOT --------------------------------------------
# Two functions, deliberately separate -- see the note above.


def extract_circles(m):
    """Pyomo model -> plain Python. After this, nothing needs the model again.

    Arguments:
        m: Pyomo concrete model, solved or merely initialized

    Returns:
        dict with the box dimensions `A` and `B`, the `radii`, and the `centers`
    """
    return {
        "A": pyo.value(m.box_height),
        "B": pyo.value(m.box_width),
        "radii": {i: pyo.value(m.R[i]) for i in m.CIRCLES},
        "centers": {i: (pyo.value(m.x[i]), pyo.value(m.y[i])) for i in m.CIRCLES},
    }


# The house hatch sequence, copied from figures/plots/_house.py as a LITERAL:
# that module is not on disk on Colab, so it must not be imported here. The
# first three entries are three different texture FAMILIES -- lines, dots,
# crosses -- so every pair of circles differs in kind rather than in the sign
# of a slope, which is the cue that survives greyscale printing.
HATCHES = ("///", "...", "xxx", "|||", "---", "\\\\\\")
SHADE_ALPHA = 0.18
OKABE_BLUE = "#0072B2"  # Okabe-Ito; the only color in the figure


def plot_circles(results, start=None, title=None, ax=None, span=None):
    """Draw one arrangement of circles in a box, in the course house style.

    Arguments:
        results: dict from `extract_circles` -- the arrangement to draw solid
        start: optionally, another such dict, drawn dotted and unfilled
            underneath to show where the solver started
        title: text above the panel
        ax: draw into an existing axes instead of creating a figure
        span: force the axis range instead of computing it from the content.
            Pass the same value to several panels when their sizes are meant
            to be compared by eye -- otherwise each panel is scaled to fit and
            two boxes of different size look identical.

    Returns:
        The matplotlib Axes.
    """
    A, B = results["A"], results["B"]
    radii, centers = results["radii"], results["centers"]

    if ax is None:
        _, ax = plt.subplots(figsize=(6, 6))
    hatch_colour = plt.rcParams["hatch.color"]

    # Axis limits from the CONTENT, not from the box: at a loose initial point
    # a circle can stick out well past the walls, and clipping it hides the
    # very thing the picture is being drawn to show.
    extent = []
    for d in (results, start):
        if d is not None:
            extent += [d["A"], d["B"]]
            for i in d["radii"]:
                cx, cy = d["centers"][i]
                extent += [cx + d["radii"][i], cy + d["radii"][i]]
    span = max(extent) if span is None else span
    pad = 0.14 * span

    # Where the solver started: dotted, unfilled, underneath everything else.
    if start is not None:
        ax.add_patch(
            mpatches.Rectangle(
                (0, 0), start["B"], start["A"], facecolor="none",
                edgecolor="0.72", linestyle=":", linewidth=1.5, zorder=1,
            )
        )
        for i in start["radii"]:
            cx, cy = start["centers"][i]
            ax.add_patch(
                mpatches.Circle(
                    (cx, cy), start["radii"][i], facecolor="none",
                    edgecolor="0.72", linestyle=":", linewidth=1.5, zorder=1,
                )
            )

    # The circles: a pale hatch that survives greyscale, a blue outline, and
    # the radius written in the middle. Circle i is labelled R_i, exactly as
    # the model indexes it.
    for k, i in enumerate(radii):
        cx, cy = centers[i]
        ax.add_patch(
            mpatches.Circle(
                (cx, cy), radii[i], facecolor="0.55", alpha=SHADE_ALPHA,
                hatch=HATCHES[k % len(HATCHES)], edgecolor=hatch_colour,
                linewidth=0.0, zorder=2,
            )
        )
        ax.add_patch(
            mpatches.Circle(
                (cx, cy), radii[i], facecolor="none", edgecolor=OKABE_BLUE,
                linewidth=2.4, zorder=3,
            )
        )
        # The label sits half a radius ABOVE the center, not on it. At a loose
        # initial point two circles can be nearly concentric -- the circles
        # overlap, because an initial point need not be feasible -- and two
        # labels on the same spot are unreadable. Offsetting by a fraction of
        # each circle's OWN radius separates them, because the radii differ.
        # The white box keeps the text legible over the hatching.
        ax.annotate(
            f"$R_{{{i}}}$", xy=(cx, cy + 0.5 * radii[i]), fontsize=13,
            ha="center", va="center", zorder=6,
            bbox=dict(facecolor="white", edgecolor="none", alpha=0.75, pad=1.0),
        )

    # The box, with its two dimensions written on it: A and B ARE the objective.
    ax.add_patch(
        mpatches.Rectangle((0, 0), B, A, facecolor="none", edgecolor="black",
                           linewidth=2.6, zorder=4)
    )
    off = -0.5 * pad
    ax.annotate("", xy=(off, 0), xytext=(off, A),
                arrowprops=dict(arrowstyle="<->", color="black", lw=1.3))
    ax.annotate(f"$A = {A:.2f}$", xy=(off, 0.5 * A), fontsize=12.5, ha="center",
                va="center", rotation=90,
                bbox=dict(facecolor="white", edgecolor="none", pad=1.0))
    ax.annotate("", xy=(0, off), xytext=(B, off),
                arrowprops=dict(arrowstyle="<->", color="black", lw=1.3))
    ax.annotate(f"$B = {B:.2f}$", xy=(0.5 * B, off), fontsize=12.5, ha="center",
                va="center", bbox=dict(facecolor="white", edgecolor="none", pad=1.0))

    ax.set_xlim(-pad, span + 0.5 * pad)
    ax.set_ylim(-pad, span + 0.5 * pad)
    ax.set_aspect("equal", adjustable="box")
    ax.set_xlabel("$x$ [m]")
    ax.set_ylabel("$y$ [m]")
    if title:
        ax.set_title(title, fontsize=13.5)
    if start is not None:
        ax.annotate("dotted: where the solver started", xy=(0.02, 0.98),
                    xycoords="axes fraction", fontsize=11, ha="left", va="top",
                    color="0.40")
    return ax


# Keep the initial point: we want to draw it again underneath the solution.
first_start = extract_circles(model)
plot_circles(first_start, title="Initial point")
plt.show()

# Specify the solver
solver = pyo.SolverFactory("ipopt")

# Solve the model
results = solver.solve(model, tee=True)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

# Print variable values
print("Name\tValue")
for c in model.component_data_objects(pyo.Var):
    print(c.name, "\t", pyo.value(c))

# Plot the solution, with the initial point dotted underneath
plot_circles(
    extract_circles(model), start=first_start, title="Solution from the first start"
)
plt.show()

# Print constraints
for c in model.component_data_objects(pyo.Constraint):
    print(
        c.name,
        "\t",
        pyo.value(c.lower),
        "\t",
        pyo.value(c.body),
        "\t",
        pyo.value(c.upper),
    )

# Initialize and print the model

# Add your solution here

# Plot initial point
# Add your solution here

# Solve the model
# Add your solution here

# Plot solution
# Add your solution here

# ---------- SOLVE ----------------------------------------------------------
# Nothing is plotted in this cell. This is the only cell below that needs a
# solver.
#
# THE INSTANCE. Biegler prints no numerical radii -- his Figure 4.4 is a sketch
# -- so R = (1.0, 0.7, 0.6) is ours. It was chosen because it has *two*
# certified local minima that are 13.6% apart in perimeter: far enough that the
# two boxes differ visibly on the page, and the radii are visibly unequal, as
# "three circles of different sizes" requires. Other triples give only one, and
# then the picture below is impossible to draw.
lecture_radii = {"1": 1.0, "2": 0.7, "3": 0.6}

# The loose box every run starts from. Biegler's second bullet on p. 68 is
# "change the initial arrangement of the circles and manually shrink the walls
# of the box", so the starting walls are part of the story and get drawn.
BOX_INIT = 4.6

# Two starting arrangements, both feasible, both describable in words.
starts = {
    "strung along the diagonal": {"1": (1.0, 1.0), "2": (2.4, 2.3), "3": (3.4, 3.4)},
    "one circle per corner": {"1": (1.2, 1.2), "2": (3.8, 0.7), "3": (0.7, 3.9)},
}

# ⚠ R_2 SITS SLIGHTLY OFF THE DIAGONAL ON PURPOSE, and that is a lesson rather
# than bookkeeping. This model is symmetric under swapping the two axes
# (A <-> B and x <-> y together). A start with all three centers exactly on the
# line y = x is a fixed point of that symmetry, and Ipopt walks straight down
# it to the symmetric stationary point: a *square* box, perimeter 14.885, which
# is not a local minimum at all. Measured with Pyomo 6.10.1 + Ipopt: centers
# (1.0, 1.0), (2.3, 2.3), (3.4, 3.4) in a box of 4.6 return 14.885, and the
# same start with the box written as 4.599999999999997 returns 11.488. The
# tie-break is being decided by the last bits of the mantissa. Break the
# symmetry yourself; do not leave it to round-off.


def solve_circle_packing(radii, centers, box_init=BOX_INIT):
    """Build, initialize from a named arrangement, and solve. Returns the model."""
    m = create_circle_model(radii)
    m.box_height = box_init
    m.box_width = box_init
    for i, (x0, y0) in centers.items():
        m.x[i] = x0
        m.y[i] = y0
    status = pyo.SolverFactory("ipopt").solve(m)
    assert pyo.check_optimal_termination(status), (
        f"solve failed: {status.solver.termination_condition}"
    )
    return m


def extract_point(m):
    """Pyomo model -> plain Python. After this, nothing needs the model again."""
    return {
        "A": pyo.value(m.box_height),
        "B": pyo.value(m.box_width),
        "perimeter": pyo.value(m.obj),
        "centers": {i: (pyo.value(m.x[i]), pyo.value(m.y[i])) for i in m.CIRCLES},
    }


points = {
    name: extract_point(solve_circle_packing(lecture_radii, centers))
    for name, centers in starts.items()
}

for name, p in points.items():
    print(
        f"start {name:28s} -> A = {p['A']:.4f} m, B = {p['B']:.4f} m, "
        f"perimeter = {p['perimeter']:.4f} m"
    )

names = list(starts)
penalty = 100.0 * (points[names[1]]["perimeter"] / points[names[0]]["perimeter"] - 1.0)
print(f"the two answers differ by {penalty:.1f}% of the better one")

# ---------- RESTRICTED MULTISTART DIAGNOSTIC -------------------------------
# Search the feasible set intersected with a ball about each returned point.
# Finite local solves do not prove that the ball contains no better point.
# A better point in this fixed ball also need not refute local minimality in
# a smaller neighborhood. The exact local certificates are derived below.

DELTA = 0.5
N_TRY = 40


def best_in_ball(radii, z_star, delta=DELTA, n_try=N_TRY, seed=17):
    """Best value found by the finite restricted multistart experiment.

    This is a search diagnostic, not a certified minimum over the ball.
    """
    rng = random.Random(seed)  # seeded: record the sampler for a reproducible diagnostic
    circles = list(radii)
    n = 2 + 2 * len(circles)
    m_units = pyo.units.m
    best = z_star["perimeter"]
    opt = pyo.SolverFactory("ipopt")

    for _ in range(n_try):
        # A point uniform in the ball: random direction, radius ~ U(0,1)^(1/n).
        d = [rng.gauss(0.0, 1.0) for _ in range(n)]
        scale = delta * rng.random() ** (1.0 / n) / sum(v * v for v in d) ** 0.5
        d = [v * scale for v in d]

        m = create_circle_model(radii)
        m.box_height = z_star["A"] + d[0]
        m.box_width = z_star["B"] + d[1]
        for k, i in enumerate(circles):
            m.x[i] = z_star["centers"][i][0] + d[2 + 2 * k]
            m.y[i] = z_star["centers"][i][1] + d[3 + 2 * k]

        # The ball itself, as one extra constraint on the same model.
        m.trust_region = pyo.Constraint(
            expr=(m.box_height - z_star["A"] * m_units) ** 2
            + (m.box_width - z_star["B"] * m_units) ** 2
            + sum(
                (m.x[i] - z_star["centers"][i][0] * m_units) ** 2
                + (m.y[i] - z_star["centers"][i][1] * m_units) ** 2
                for i in circles
            )
            <= delta**2 * m_units**2
        )

        status = opt.solve(m)
        if pyo.check_optimal_termination(status):
            best = min(best, pyo.value(m.obj))
    return best


ball_check = {}
for name, p in points.items():
    ball_check[name] = best_in_ball(lecture_radii, p)
    verdict = "NO LOWER VALUE FOUND" if ball_check[name] > p["perimeter"] - 1e-3 else "LOWER VALUE FOUND IN TEST BALL"
    print(
        f"{name:28s} perimeter {p['perimeter']:.5f}, best in ball of radius "
        f"{DELTA} = {ball_check[name]:.5f}  ->  {verdict}"
    )

# ---------- ARCHIVE --------------------------------------------------------
# figures/results/packing-local-solutions.json, committed to the repository.
# `source_tag` points at the model cell above, so scripts/check_results_fresh.py
# can tell you when the model has changed and these numbers have not.
# A no-op on Colab, where there is nothing to commit to.

results = {
    "radii": dict(lecture_radii),
    "box_init": BOX_INIT,
    "delta": DELTA,
    "panels": [
        {
            "name": name,
            "start": {i: list(starts[name][i]) for i in starts[name]},
            "A": points[name]["A"],
            "B": points[name]["B"],
            "perimeter": points[name]["perimeter"],
            "centers": {i: list(points[name]["centers"][i]) for i in lecture_radii},
            "best_in_ball": ball_check[name],
        }
        for name in starts
    ],
}

helper.save_results(
    "packing-local-solutions",
    results,
    notebook="notebooks/1-dev/NLP.ipynb",
    source_tag="handout:circle-packing-model",
    source_cells=[("notebooks/1-dev/NLP.ipynb", t) for t in ['handout:circle-packing-model', 'results:packing-pair-solve', 'results:packing-pair-diagnostic', 'results:packing-pair-extract']],
    description="Biegler Example 4.4 with R = (1.0, 0.7, 0.6), solved twice "
    "from two starting arrangements: two returned arrangements 13.6% apart in "
    "perimeter, with a finite restricted multistart check of radius 0.5.",
    solver="Ipopt via Pyomo",
)

# The PLOTTING function: it takes the extracted results, not the Pyomo model.
#
# Everything it needs is in `results`, and nothing else -- no model, no solver.
# That is what lets you re-run this one cell as many times as it takes to get
# the labels where you want them.
#
# This cell is tagged `figure:packing-local-solutions`, which makes it the
# single source of the figure the lecture handout prints. See figures/README.md;
# figures/render_from_notebook.py re-runs exactly this cell against the archived
# JSON when the house style changes. That driver binds `helper` for the cell,
# so this cell may use it without importing anything.

from matplotlib.patches import Circle, Rectangle


def plot_packing_local_solutions(results):
    """Two local solutions of the circle-packing NLP, from two starting guesses.

    WHY THIS FIGURE EXISTS. The model above is nonconvex -- the No Overlap
    constraints are -- so the convexity theorem does not apply and

        "the optimum you find depends on where you start."

    This is the evidence for that sentence, and it is Biegler's own second
    bullet on p. 68. Same model, same solver, same tolerances, two starting
    arrangements, two different answers.

    The restricted multistart experiment is numerical evidence, not a proof.
    The separate exact active-equation/multiplier derivation below gives local
    certificates for the idealized versions of the plotted arrangements.
    Global optimality is a separate question.

    GREYSCALE. The three circles carry three genuinely different hatch families
    -- lines, dots, crosses -- and their own labels, so they are told apart
    without color. The starting arrangement is dotted and unfilled against the
    solid, filled solution: a linestyle contrast, not a color one. The box is
    black in both panels.
    """
    radii = results["radii"]
    panels = results["panels"]
    box0 = results["box_init"]
    circles = list(radii)

    # Hatch literals copied from figures/plots/_house.py, which is NOT on disk
    # on Colab and must not be imported here:
    #   HATCH_CYCLE = ("///", "\\\\\\", "...", "xxx", "|||", "---")
    #   SHADE_ALPHA = 0.18
    # ⚠ HATCH_CYCLE[1], "\\\", is skipped deliberately. Three circles need
    # three textures, and these three spend three different texture FAMILIES,
    # so every pair differs in kind. Restoring [1] puts two of the three in the
    # same family, separated only by the SIGN of the slope -- the most fragile
    # cue here, because the hatch is drawn pale at low alpha inside circles
    # that touch.
    hatches = ("///", "...", "xxx")
    shade_alpha = 0.18
    hatch_colour = plt.rcParams["hatch.color"]
    blue = "#0072B2"  # Okabe-Ito; the solved circles are the only color here
    lim_lo, lim_hi = -0.62, 5.00

    def draw_circles(ax, centers, solution):
        """One arrangement. `solution` toggles solved vs starting guess."""
        for k, i in enumerate(circles):
            x, y = centers[i]
            if solution:
                ax.add_patch(Circle((x, y), radii[i], facecolor="0.55",
                                    alpha=shade_alpha, hatch=hatches[k],
                                    edgecolor=hatch_colour, linewidth=0.0,
                                    zorder=2))
                ax.add_patch(Circle((x, y), radii[i], facecolor="none",
                                    edgecolor=blue, linewidth=2.4, zorder=3))
                ax.annotate("$R_%d$" % (k + 1), xy=(x, y), fontsize=13,
                            ha="center", va="center", zorder=6)
            else:
                ax.add_patch(Circle((x, y), radii[i], facecolor="none",
                                    edgecolor="0.72", linestyle=":",
                                    linewidth=1.5, zorder=1))

    def draw_box(ax, A, B):
        """The enclosing box, with its two dimensions written on it."""
        ax.add_patch(Rectangle((0.0, 0.0), B, A, facecolor="none",
                               edgecolor="black", linewidth=2.6, zorder=4))
        ax.annotate("", xy=(-0.30, 0.0), xytext=(-0.30, A),
                    arrowprops=dict(arrowstyle="<->", color="black", lw=1.3))
        ax.annotate("$A = %.2f$" % A, xy=(-0.30, 0.5 * A), fontsize=12.5,
                    ha="center", va="center", rotation=90,
                    bbox=dict(facecolor="white", edgecolor="none", pad=1.0))
        ax.annotate("", xy=(0.0, -0.30), xytext=(B, -0.30),
                    arrowprops=dict(arrowstyle="<->", color="black", lw=1.3))
        ax.annotate("$B = %.2f$" % B, xy=(0.5 * B, -0.30), fontsize=12.5,
                    ha="center", va="center",
                    bbox=dict(facecolor="white", edgecolor="none", pad=1.0))

    def frame(ax, title):
        """Apply common axes limits, labels, and title to one panel."""
        ax.set_xlim(lim_lo, lim_hi)
        ax.set_ylim(lim_lo, lim_hi)
        ax.set_aspect("equal", adjustable="box")
        ax.set_xlabel("$x$")
        ax.set_title(title, fontsize=13.5)
        ax.set_xticks([0, 1, 2, 3, 4])
        ax.set_yticks([0, 1, 2, 3, 4])

    # ⚠ figsize is set HERE. helper.set_plotting_style() overrides the house
    # default for on-screen readability and render_from_notebook.py does not,
    # so a figure that relied on the default would come out at two different
    # aspect ratios depending on who generated it.
    fig, axes = plt.subplots(1, 2, figsize=(9.6, 5.0))
    axes[0].set_ylabel("$y$")

    penalty = 100.0 * (panels[1]["perimeter"] / panels[0]["perimeter"] - 1.0)
    # ⚠ f-strings, not %-formatting, for these two: scripts/check_greyscale.py
    # blanks any line whose first non-blank character is `%` (it reads it as a
    # Jupyter line magic), so a `% (...)` continuation line makes the whole
    # notebook score "unparseable".
    titles = [
        f"start: {panels[0]['name']}\n$2(A+B) = {panels[0]['perimeter']:.2f}$",
        f"start: {panels[1]['name']}\n$2(A+B) = {panels[1]['perimeter']:.2f}$"
        f"   $(+{penalty:.1f}\\%)$",
    ]

    for ax, panel, title in zip(axes, panels, titles):
        frame(ax, title)
        # The loose box every run starts from: p. 68 is about shrinking these
        # walls, so the walls have to be on the page.
        ax.add_patch(Rectangle((0.0, 0.0), box0, box0, facecolor="none",
                               edgecolor="0.72", linestyle=":", linewidth=1.5,
                               zorder=1))
        draw_circles(ax, panel["start"], solution=False)
        draw_circles(ax, panel["centers"], solution=True)
        draw_box(ax, panel["A"], panel["B"])

    axes[0].annotate("dotted: where the solver started", xy=(0.10, 4.30),
                     fontsize=11.5, ha="left", va="top", color="0.40")

    fig.tight_layout()
    for axis in fig.axes:
        axis.tick_params(colors="black")
        for spine in axis.spines.values():
            spine.set_color("black")
    return fig


fig = plot_packing_local_solutions(results)

# Write media/figures/packing-local-solutions.{png,pdf} -- what the lecture
# handout \includegraphics. A no-op on Colab, where there is no repo to write to.
helper.save_figure(fig, "packing-local-solutions")

def multistart_circle_packing(
    circle_data, n_starts=50, seed=0, spread=40, box_init=40, tee=False
):
    """Solve the circle packing model from many random initial points

    Arguments:
        circle_data: dictionary with keys=circle name and value=radius (float)
        n_starts: number of random initializations (int)
        seed: seed for the random number generator, so the result reproduces
        spread: width of the uniform window the centers are drawn from
        box_init: initial value for both box dimensions
        tee: print solver output

    Returns:
        results: pandas DataFrame, one row per start, with the objective and
            the box dimensions (NaN if that start did not solve)
        arrangements: dict keyed by start, each holding the extracted `start`
            point and the extracted `solution` (None if that start failed).
            This is what lets us go back and LOOK at a chosen solution; the
            DataFrame alone records that the local solutions differ, not what
            any of them looks like.
    """
    # Reproducibility. Reseeding here means this cell gives the same answer no
    # matter which other cells were run first.
    random.seed(seed)

    solver = pyo.SolverFactory("ipopt")

    rows = []
    arrangements = {}
    for k in range(n_starts):
        m = create_circle_model(circle_data)
        initialize_circle_model(
            m, height_init=box_init, width_init=box_init, spread=spread
        )
        # Keep the initial point BEFORE solving -- afterwards it is gone, and
        # for a nonconvex model the initial point is half of what makes an
        # answer reportable.
        # ⚠ Extracting it draws no random numbers, so the sampling sequence,
        # and therefore every number this section prints, is unchanged.
        start_point = extract_circles(m)

        res = solver.solve(m, tee=tee)
        if pyo.check_optimal_termination(res):
            arrangements[k] = {"start": start_point, "solution": extract_circles(m)}
            rows.append(
                {
                    "start": k,
                    "perimeter": pyo.value(m.obj),
                    "A": pyo.value(m.box_height),
                    "B": pyo.value(m.box_width),
                }
            )
        else:
            arrangements[k] = {"start": start_point, "solution": None}
            rows.append(
                {"start": k, "perimeter": float("nan"), "A": float("nan"),
                 "B": float("nan")}
            )

    return pd.DataFrame(rows), arrangements


# The multistart experiment, named ONCE. The archive cell below records these
# four numbers alongside the perimeters, and a literal repeated in two cells is
# a literal that will eventually disagree with itself.
MS_STARTS = 50
MS_SEED = 0
MS_SPREAD = 40
MS_BOX_INIT = 40

ms, ms_points = multistart_circle_packing(
    circle_data,
    n_starts=MS_STARTS,
    seed=MS_SEED,
    spread=MS_SPREAD,
    box_init=MS_BOX_INIT,
)

n_failed = int(ms["perimeter"].isna().sum())
print(f"{len(ms)} starts, {n_failed} failed to converge")
print(f"Best perimeter found: {ms['perimeter'].min():.4f}")

# How many distinct local solutions did we land on?
counts = (
    ms["perimeter"]
    .round(3)
    .value_counts()
    .rename_axis("perimeter")
    .reset_index(name="number of starts")
    .sort_values("perimeter")
    .reset_index(drop=True)
)
counts["fraction"] = counts["number of starts"] / len(ms)
display(counts)

# ---------- EXTRACT + ARCHIVE ---------------------------------------------
# pandas -> plain Python, then figures/results/packing-multistart-histogram.json,
# committed to the repository. After this cell nothing below touches `ms`, a
# Pyomo model, or a solver: the plotting cell reads `results` and nothing else,
# so the histogram can be redrawn as many times as it takes without paying for
# fifty solves again.
#
# `source_tag` points at the model cell above, so scripts/check_results_fresh.py
# can tell you when the model has changed and these numbers have not. Writing is
# a no-op on Colab, where there is nothing to commit to.
#
# ⚠ WHAT IS ARCHIVED IS THE FIFTY PERIMETERS, NOT THE HISTOGRAM'S BIN HEIGHTS.
# The number of bins is a plotting decision and belongs in the plot cell;
# freezing the bins here would make a *style* change to this figure impossible
# without re-solving, and avoiding exactly that is the archive's only job.
# A start that failed to converge is stored as `null` rather than dropped, so
# the archive still records how many starts were run.

results = {
    "radii": dict(circle_data),
    "n_starts": MS_STARTS,
    "seed": MS_SEED,
    "spread": MS_SPREAD,
    "box_init": MS_BOX_INIT,
    "perimeters": [None if pd.isna(v) else float(v) for v in ms["perimeter"]],
    "best": float(ms["perimeter"].min()),
    # The `counts` table printed above: one row per distinct local solution.
    "spikes": helper.table(counts),
}

helper.save_results(
    "packing-multistart-histogram",
    results,
    notebook="notebooks/1-dev/NLP.ipynb",
    source_tag="handout:circle-packing-model",
    source_cells=[("notebooks/1-dev/NLP.ipynb", t) for t in ['handout:circle-packing-model', 'results:packing-initializer', 'results:packing-notebook-data', 'results:packing-multistart-solve', 'results:packing-multistart-count', 'results:packing-multistart-extract']],
    description=(
        "Circle packing on the notebook's instance R = (10, 5, 3): perimeters "
        "from 50 random starts, seed 0, centers drawn uniformly from a window "
        "of width 40. The perimeters fall in three returned objective groups "
        "rather than a spread."
    ),
    solver="Ipopt via Pyomo",
)

# The PLOTTING function: it takes the extracted results, not the Pyomo model.
#
# This cell is tagged `figure:packing-multistart-histogram`, which makes it the
# single source of the histogram the lecture handout prints. See
# figures/README.md; figures/render_from_notebook.py re-runs exactly this cell
# against the archived JSON when the house style changes, so there is no second
# copy of the plotting code anywhere. That driver binds `helper` for the cell,
# so this cell may use it without importing anything.


def plot_packing_multistart_histogram(results):
    """Perimeters from many random starts: a few sharp spikes, not a spread.

    WHY THIS FIGURE EXISTS. The two-panel figure earlier in this notebook shows
    that *two* starting arrangements give two different boxes. It cannot show
    how many different boxes there are, or how often the solver finds the best
    one. This does: every start is one sample, and the samples pile up on a
    handful of values because Ipopt converges to a small number of distinct
    local solutions.

    ⚠ THIS IS A LARGER INSTANCE THAN THE TWO-PANEL FIGURE, on purpose --
    R = (10, 5, 3) here against R = (1.0, 0.7, 0.6) there -- so the perimeters
    are around 98 rather than around 11. The two figures are two different
    experiments and must not be read as one; the axis label and the title both
    carry the instance so a reader cannot mistake them.

    THE COUNTS ARE PRINTED, and that is reproducible only because the seed is.
    An unseeded sampler gives a different count every run; with the seed quoted,
    "42 of 50" is a number someone else can obtain. The *best value* is far more
    robust than the count -- it is the same for every seed tried, while the
    split between the spikes is not -- so the best value is drawn as a rule
    across the whole axes and the counts are annotations on the bars.

    GREYSCALE. There is one data series and it is achromatic: grey bars with
    black edges. The only color is the vermillion rule at the best value
    found, which is also DASHED, so it is identified without color.
    """
    perimeters = [p for p in results["perimeters"] if p is not None]
    n_starts = results["n_starts"]
    seed = results["seed"]
    best = results["best"]
    radii = results["radii"]

    # ⚠ figsize is set HERE. helper.set_plotting_style() overrides the house
    # default for on-screen readability and render_from_notebook.py does not,
    # so a figure that relied on the default would come out at two different
    # aspect ratios depending on who generated it.
    fig, ax = plt.subplots(figsize=(7.4, 4.4))

    heights, edges, _ = ax.hist(perimeters, bins=30, color="0.60",
                                edgecolor="black", linewidth=1.0, zorder=2)

    # Okabe-Ito vermillion, from figures/dowling.mplstyle. Plain "red" is not in
    # the palette: it is the color a red-green color-blind reader cannot
    # separate from the black bar edges, and it prints as almost the same grey.
    ax.axvline(best, color="#D55E00", linestyle="--", linewidth=2.4, zorder=3,
               label=f"best found = {best:.3f} m")

    # One annotation per spike, read off the archived counts table rather than
    # typed in: the number of spikes depends on the seed and the sampler, so a
    # hard-coded three would break by drawing the wrong labels rather than by
    # raising an error.
    cols = results["spikes"]["columns"]
    ip, ic = cols.index("perimeter"), cols.index("number of starts")
    rows = results["spikes"]["rows"]

    # ⚠ EACH LABEL IS PLACED OVER ITS BAR, NOT OVER ITS PERIMETER. A spike sits
    # somewhere inside a bin, and a bin here is about 0.26 m wide, so a label
    # centered on the perimeter itself hangs off the side of the bar it belongs
    # to -- and for the leftmost spike, which is the whole point of the figure,
    # off the left edge of the axes entirely. Found by rendering and looking.
    for row in rows:
        k = min(range(len(heights)), key=lambda j: abs(row[ip] - edges[j]))
        k = min(k, len(heights) - 1)
        if row[ip] < edges[k]:
            k = max(k - 1, 0)
        ax.annotate(
            f"{int(row[ic])} of {n_starts}",
            xy=(0.5 * (edges[k] + edges[k + 1]), heights[k]),
            xytext=(0, 5),
            textcoords="offset points",
            ha="center",
            va="bottom",
            fontsize=11.5,
            zorder=4,
            # The label over the tallest bar is crossed by the dashed rule at
            # the best value; a white patch behind the text keeps both legible
            # without moving either off the thing it labels.
            bbox=dict(facecolor="white", edgecolor="none", pad=0.8),
        )

    # Headroom above the tallest bar for its label, and a margin at each end so
    # the outer two labels are not clipped by the frame.
    ax.set_ylim(0, 1.17 * max(heights))
    pad = 0.05 * (edges[-1] - edges[0])
    ax.set_xlim(edges[0] - pad, edges[-1] + pad)

    r_text = ", ".join(f"{radii[i]:g}" for i in radii)
    ax.set_xlabel("Perimeter of the enclosing box, $2(A + B)$ [m]")
    ax.set_ylabel("Number of starts")
    ax.set_title(
        f"{n_starts} random starts, seed {seed};  $R = ({r_text})$", fontsize=13.5
    )
    ax.legend(loc="upper right")

    fig.tight_layout()
    for axis in fig.axes:
        axis.tick_params(colors="black")
        for spine in axis.spines.values():
            spine.set_color("black")
    return fig


fig = plot_packing_multistart_histogram(results)

# Write media/figures/packing-multistart-histogram.{png,pdf} -- what the lecture
# handout \includegraphics. A no-op on Colab, where there is no repo to write to.
helper.save_figure(fig, "packing-multistart-histogram")

# Which start gave the best box?
best_start = int(ms["perimeter"].idxmin())
print(f"Best start: #{best_start}")
print(f"  perimeter = {ms.loc[best_start, 'perimeter']:.4f} m")
print(f"  A = {ms.loc[best_start, 'A']:.4f} m, B = {ms.loc[best_start, 'B']:.4f} m")

# One representative start per spike, chosen programmatically from `counts`.
spikes = ms["perimeter"].round(3)
representatives = [
    {
        "perimeter": p,
        "n_starts": int(n),
        # The first start that landed on this spike. Any member would do; the
        # first is the reproducible choice.
        "start": int(spikes[spikes == p].index[0]),
    }
    for p, n in zip(counts["perimeter"], counts["number of starts"])
]

for r in representatives:
    print(
        f"perimeter {r['perimeter']:.3f} m  <-  {r['n_starts']:2d} of {len(ms)} "
        f"starts, e.g. start #{r['start']}"
    )

# One shared scale across the panels, so the boxes can be compared by eye.
# Without this each panel is fitted to its own contents and the biggest box
# looks the same size as the smallest.
shared_span = max(
    max(ms_points[r["start"]]["solution"]["A"], ms_points[r["start"]]["solution"]["B"])
    for r in representatives
)

fig, axes = plt.subplots(
    1, len(representatives), figsize=(4.6 * len(representatives), 5.2), squeeze=False
)
for ax, r in zip(axes[0], representatives):
    best = " (best found)" if r is representatives[0] else ""
    plot_circles(
        ms_points[r["start"]]["solution"],
        title=(
            f"{r['n_starts']} of {len(ms)} starts{best}\n"
            f"$2(A+B) = {r['perimeter']:.3f}$ m"
        ),
        ax=ax,
        span=shared_span,
    )
    if ax is not axes[0][0]:
        ax.set_ylabel("")
fig.tight_layout()
plt.show()