Matplotlib Optimization Visualization

Matplotlib Visualization for Optimization

You are an expert in scientific visualization for operations research. This skill covers the standard figure types of computational optimization papers — convergence curves, Gantt charts, route plots, Pareto fronts, and performance profiles — plus the publication-quality mechanics (vector output, font matching, column-width sizing) that journals and conferences require. Use the pattern catalog below: each pattern gives the motivation, a complete implementation, and the pitfall that most often ruins the figure.

Initial Assessment

Before producing any figure, establish:

• Venue and column geometry. Single-column (~3.3–3.5 in) or double-column (~7 in) figure? IEEE, INFORMS, Springer, and Elsevier column widths differ; the figure must be designed at its final printed width.
• Target format. PDF or EPS vector for the paper; PNG only for previews, slides, or raster-heavy panels. Some journals still require EPS or TIFF — confirm before styling.
• Data shape. Is the experiment data already in tidy form (one row per run, columns for instance, algorithm, seed, time, objective)? If not, fix the table first — plotting code should never reshape ad hoc.
• What the figure must argue. Anytime behavior (convergence plot), final quality distribution (box/strip plot), robustness across instances (performance profile), structural correctness (route/Gantt plot), or trade-offs (Pareto plot). One claim per figure.
• Minimization or maximization. Determines np.minimum.accumulate vs np.maximum.accumulate, axis direction, and which corner of a Pareto plot is "good."
• Number of seeds and instances. Fewer than ~5 seeds: plot individual runs, not a band. Many instances: aggregate with performance profiles, not 50 separate convergence plots.
• Time axis semantics. Wall-clock seconds, CPU seconds, or evaluation count? Mixed hardware makes wall-clock comparisons unfair; evaluation counts hide per-iteration cost differences. State the choice on the axis label.
• Scale needs. Objectives spanning orders of magnitude, or late-stage differences of None:

"""Set global rcParams for camera-ready figures. Call once, before creating figures.""" mpl.rcParams.update({ "font.family": "serif", "font.size": basept, "axes.labelsize": basept, "axes.titlesize": basept, "xtick.labelsize": basept - 1, "ytick.labelsize": basept - 1, "legend.fontsize": basept - 1, "mathtext.fontset": "cm", # Computer Modern: matches LaTeX math "axes.spines.top": False, "axes.spines.right": False, "axes.grid": True, "grid.linewidth": 0.4, "grid.alpha": 0.4, "lines.linewidth": 1.2, "lines.markersize": 3.5, "legend.frameon": False, "pdf.fonttype": 42, # embed TrueType outlines, never Type 3 "ps.fonttype": 42, "savefig.dpi": 300, # raster fallback resolution "figure.constrained_layout.use": True, })

def figsize(widthin: float, aspect: float = 0.62) -> tuple[float, float]: """Figure size at FINAL printed width; aspect = height/width (0.62 ~ golden ratio).""" return (widthin, round(widthin * aspect, 2))

setpublicationstyle() print(figsize(JOURNALWIDTHSIN["singlecolumn"]))

Expected: (3.5, 2.17) -- include at natural size in LaTeX so 8 pt text prints at 8 pt

**Pitfall:** Setting `figsize=(10, 6)` and later scaling with `width=\columnwidth` in LaTeX shrinks the figure by ~65%, turning 10 pt fonts into ~3.5 pt fonts. This is the single most common defect in submitted optimization papers. Fix the width in inches here and never rescale downstream.

### Pattern 2 — Save vector master plus raster preview

The paper needs a vector PDF (crisp at any zoom, searchable text); day-to-day inspection needs a PNG that opens fast and embeds in notebooks and chat. Save both from one call, and close the figure so long experiment scripts do not leak memory through hundreds of open canvases.

```python
from pathlib import Path

import matplotlib
matplotlib.use("Agg")          # headless backend: works on servers and in CI
import matplotlib.pyplot as plt
from matplotlib.figure import Figure

def save_figure(fig: Figure, stem: str | Path,
                formats: tuple[str, ...] = ("pdf", "png")) -> list[Path]:
    """Save one figure in several formats; pdf is the paper artifact, png the preview."""
    stem = Path(stem)
    stem.parent.mkdir(parents=True, exist_ok=True)
    written: list[Path] = []
    for ext in formats:
        target = stem.with_suffix(f".{ext}")
        fig.savefig(target)    # constrained_layout already handles spacing
        written.append(target)
    plt.close(fig)
    return written

fig, ax = plt.subplots(figsize=(3.5, 2.17))
ax.plot([0, 1, 2], [3, 1, 2])
ax.set_xlabel("iteration")
ax.set_ylabel("objective")
paths = save_figure(fig, "figures/demo_curve")
print([p.name for p in paths])
# Expected: ['demo_curve.pdf', 'demo_curve.png'] -- vector master plus raster preview

Pitfall: bbox_inches="tight" recomputes the canvas size at save time, so the saved figure is not the width you designed, and a grid of "3.5 in" figures ends up with three slightly different widths. With constrained_layout enabled (Pattern 1) you do not need tight; if you must crop, accept that printed width changed and re-check font sizes.

Convergence Patterns

Pattern 3 — Best-so-far incumbent as a step curve

The incumbent objective is piecewise constant between improvements. Reduce the raw evaluation log to improvement events, then draw a right-continuous staircase and extend the final incumbent to the end of the run, otherwise the curve visually "stops" at the last improvement and hides the long tail without progress.

import numpy as np
import matplotlib.pyplot as plt

def best_so_far(times: np.ndarray, objectives: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Reduce a raw evaluation log to incumbent improvement events (minimization)."""
    order = np.argsort(times, kind="stable")
    t = times[order]
    f = np.minimum.accumulate(objectives[order])
    keep = np.ones(t.size, dtype=bool)
    keep[1:] = f[1:]  None:
    """Right-continuous step curve of the incumbent objective over time."""
    t, f = best_so_far(times, objectives)
    t_ext = np.append(t, run_end)      # hold the last incumbent to the end of the run
    f_ext = np.append(f, f[-1])
    ax.step(t_ext, f_ext, where="post", label=label)

rng = np.random.default_rng(7)
t = np.sort(rng.uniform(0.0, 60.0, size=400))
f = 1000.0 * np.exp(-t / 25.0) + rng.normal(0.0, 15.0, size=400) + 200.0
fig, ax = plt.subplots(figsize=(3.5, 2.17))
plot_incumbent(ax, t, f, run_end=60.0, label="ILS")
ax.set_xlabel("wall-clock time (s)")
ax.set_ylabel("best objective found")
ax.legend()
fig.savefig("incumbent_demo.pdf")
# Expected: a monotone non-increasing staircase from ~1150 down to ~210, flat after the
# last improvement and extended to t = 60

Pitfall: Logging every evaluation instead of every improvement makes log files of metaheuristics explode (millions of rows) and plotting slow. Log improvement events only — (time, new_incumbent) — plus one final row at the time limit. The plot needs nothing else.

Pattern 4 — Quantile band over seeds on a common time grid

Single-seed curves overstate differences; averages over seeds require all trajectories sampled at the same time points. Resample each seed's staircase onto a shared grid with previous-value interpolation (never linear — see Pattern 3), then plot the median with a 25–75% band. A log-spaced grid gives the early phase, where most improvement happens, enough resolution.

import numpy as np
import matplotlib.pyplot as plt

def step_resample(event_t: np.ndarray, event_f: np.ndarray,
                  grid: np.ndarray) -> np.ndarray:
    """Sample a step trajectory on a grid with previous-value (staircase) interpolation."""
    idx = np.searchsorted(event_t, grid, side="right") - 1
    out = np.full(grid.size, np.nan)   # NaN before the first incumbent exists
    seen = idx >= 0
    out[seen] = event_f[idx[seen]]
    return out

def plot_convergence_band(ax: plt.Axes, runs: list[tuple[np.ndarray, np.ndarray]],
                          grid: np.ndarray, label: str, color: str) -> None:
    """Median incumbent over seeds with a 25-75% quantile band on a shared time grid."""
    curves = np.vstack([step_resample(t, f, grid) for t, f in runs])
    q25, q50, q75 = np.nanpercentile(curves, [25.0, 50.0, 75.0], axis=0)
    ax.plot(grid, q50, color=color, label=label)
    ax.fill_between(grid, q25, q75, color=color, alpha=0.25, linewidth=0)

def synthetic_runs(rate: float, noise: float, n_seeds: int,
                   seed: int) -> list[tuple[np.ndarray, np.ndarray]]:
    """Generate improvement logs (event times, incumbent values) for one algorithm."""
    rng = np.random.default_rng(seed)
    runs: list[tuple[np.ndarray, np.ndarray]] = []
    for _ in range(n_seeds):
        n_events = int(rng.integers(20, 40))
        t = np.sort(rng.uniform(0.05, 60.0, size=n_events))
        f = 500.0 * np.exp(-rate * t) + 100.0 + rng.normal(0.0, noise, size=n_events)
        runs.append((t, np.minimum.accumulate(f)))
    return runs

grid = np.geomspace(0.1, 60.0, num=200)    # log-spaced: early progress gets resolution
fig, ax = plt.subplots(figsize=(3.5, 2.17))
plot_convergence_band(ax, synthetic_runs(0.10, 8.0, 10, seed=1), grid, "ALNS", "C0")
plot_convergence_band(ax, synthetic_runs(0.06, 8.0, 10, seed=2), grid, "GA", "C1")
ax.set_xscale("log")
ax.set_xlabel("wall-clock time (s)")
ax.set_ylabel("best objective (median, IQR over 10 seeds)")
ax.legend()
fig.savefig("convergence_band.pdf")
# Expected: ALNS band drops faster and sits below the GA band after ~5 s; bands overlap
# early, which honestly shows the early phase is not statistically separated

Pitfall: Computing the mean instead of the median lets one bad seed drag the whole curve, and computing quantiles over runs of different lengths without the NaN handling above silently mixes "no incumbent yet" with real values. Use np.nanpercentile, and never extrapolate a seed's trajectory beyond its own run end — if seeds have different budgets, truncate the grid to the shortest budget.

Solution-Structure Patterns

Pattern 5 — Route plots for TSP/VRP solutions

A route plot is the fastest sanity check for routing output: crossing edges suggest missed 2-opt moves, a giant route next to tiny ones suggests broken capacity handling. One color per vehicle, depot as a distinct marker, equal aspect so geometry is not distorted, and direction shown by an arrow on the first leg.

import numpy as np
import matplotlib.pyplot as plt

def plot_routes(ax: plt.Axes, coords: np.ndarray,
                routes: list[list[int]], depot: int = 0) -> None:
    """Draw vehicle routes over customer coordinates; one color per route, depot square."""
    cmap = plt.get_cmap("tab10")
    for k, route in enumerate(routes):
        seq = np.array([depot, *route, depot])
        xy = coords[seq]
        color = cmap(k % 10)
        ax.plot(xy[:, 0], xy[:, 1], "-", color=color, linewidth=1.0,
                label=f"route {k + 1} ({len(route)} stops)", zorder=1)
        mid = 0.5 * (xy[0] + xy[1])    # arrow on the first leg shows direction
        ax.annotate("", xy=tuple(mid), xytext=tuple(xy[0]),
                    arrowprops={"arrowstyle": "-|>", "color": color, "lw": 1.0})
    customers = np.setdiff1d(np.arange(len(coords)), [depot])
    ax.scatter(coords[customers, 0], coords[customers, 1], s=12, color="0.25", zorder=2)
    ax.scatter(coords[depot, 0], coords[depot, 1], marker="s", s=45,
               color="black", zorder=3, label="depot")
    ax.set_aspect("equal")
    ax.set_xticks([])
    ax.set_yticks([])
    ax.legend(loc="upper left", bbox_to_anchor=(1.01, 1.0))

rng = np.random.default_rng(42)
coords = rng.uniform(0.0, 100.0, size=(13, 2))
coords[0] = (50.0, 50.0)
routes = [[1, 4, 7, 10], [2, 5, 8, 11], [3, 6, 9, 12]]
fig, ax = plt.subplots(figsize=(3.5, 3.0))
plot_routes(ax, coords, routes)
fig.savefig("routes_demo.pdf")
# Expected: three colored loops through the black depot square at (50, 50), customer
# dots in gray, route legend placed outside the axes on the right

Pitfall: Forgetting set_aspect("equal") stretches the plane, so routes that cross look fine and vice versa — the one thing the figure exists to show becomes unreliable. Also resist plotting node indices for instances beyond ~30 nodes; the labels turn the plot into noise. Annotate only nodes you discuss in the text.

Pattern 6 — Gantt charts for machine schedules

A Gantt chart shows feasibility (no overlapping bars in a lane), idle time (gaps), and the makespan (dashed line) at a glance. One horizontal lane per machine, one color per job so precedence chains are traceable across machines, and bar labels only where they fit.

import matplotlib.pyplot as plt
from matplotlib.patches import Patch

Operation = tuple[int, int, float, float]      # (job, machine, start, end)

def plot_gantt(ax: plt.Axes, ops: list[Operation], n_machines: int) -> None:
    """Machine-row Gantt: one lane per machine, one color per job, makespan line."""
    cmap = plt.get_cmap("tab20")
    jobs = sorted({job for job, _, _, _ in ops})
    color_of = {job: cmap(i % 20) for i, job in enumerate(jobs)}
    for job, machine, start, end in ops:
        ax.barh(machine, end - start, left=start, height=0.6,
                color=color_of[job], edgecolor="black", linewidth=0.4)
        if end - start >= 4.0:                 # label only bars wide enough to read
            ax.text(0.5 * (start + end), machine, f"J{job}",
                    ha="center", va="center", fontsize=6)
    makespan = max(end for _, _, _, end in ops)
    ax.axvline(makespan, color="0.3", linestyle="--", linewidth=0.8)
    ax.set_yticks(range(n_machines), [f"M{m}" for m in range(n_machines)])
    ax.invert_yaxis()                          # machine 0 on top: reading order
    ax.set_xlabel("time")
    handles = [Patch(facecolor=color_of[j], edgecolor="black", label=f"job {j}")
               for j in jobs]
    ax.legend(handles=handles, ncol=min(len(jobs), 4), loc="upper center",
              bbox_to_anchor=(0.5, -0.28))

ops: list[Operation] = [
    (0, 0, 0.0, 5.0), (0, 1, 5.0, 9.0), (0, 2, 9.0, 16.0),
    (1, 1, 0.0, 5.0), (1, 0, 5.0, 11.0), (1, 2, 16.0, 21.0),
    (2, 2, 0.0, 8.0), (2, 0, 11.0, 18.0), (2, 1, 18.0, 24.0),
]
fig, ax = plt.subplots(figsize=(3.5, 2.4))
plot_gantt(ax, ops, n_machines=3)
fig.savefig("gantt_demo.pdf")
# Expected: 3-lane chart with no overlap inside any lane, dashed makespan line at
# t = 24, three-entry job legend centered below the axes

Pitfall: Coloring by machine instead of by job makes every bar in a lane the same color, which hides exactly the information a job-shop Gantt must show: how each job flows across machines. Color by job; the lane already encodes the machine. For more than ~20 jobs, drop the legend and per-bar labels, and color by a job attribute instead (due-date tightness, tardiness, family).

Comparison Patterns

Pattern 7 — Pareto front with dominated background

For bi-objective minimization, show all evaluated points in light gray and overlay the non-dominated set as a staircase. The gray cloud gives scale ("how much of the search was wasted?"); the staircase, drawn with where="post", encodes the attainment boundary exactly — straight lines between front points claim trade-off solutions that do not exist.

import numpy as np
import matplotlib.pyplot as plt

def pareto_mask(points: np.ndarray) -> np.ndarray:
    """Boolean mask of non-dominated points for bi-objective minimization."""
    order = np.lexsort((points[:, 1], points[:, 0]))   # by f1, tie-break

…

## Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

- **Author:** [hajibabaie](https://github.com/hajibabaie)
- **Source:** [hajibabaie/combinatorial-optimization-skills](https://github.com/hajibabaie/combinatorial-optimization-skills)
- **License:** MIT

Install and usage instructions live in the source repository linked above.