function load_notebook_bootstrap()
    candidates = (
        joinpath(pwd(), "scripts", "notebook_bootstrap.jl"),
        joinpath(pwd(), "..", "scripts", "notebook_bootstrap.jl"),
        joinpath(pwd(), "QuIP", "scripts", "notebook_bootstrap.jl"),
        joinpath("/content", "QuIP", "scripts", "notebook_bootstrap.jl"),
    )

    for candidate in candidates
        if isfile(candidate)
            include(candidate)
            return nothing
        end
    end

    in_colab = haskey(ENV, "COLAB_RELEASE_TAG") || haskey(ENV, "COLAB_JUPYTER_IP") || isdir(joinpath("/content", "sample_data"))
    if in_colab
        repo_dir = get(ENV, "QUIP_REPO_DIR", joinpath(pwd(), "QuIP"))
        if !isdir(repo_dir)
            println("[bootstrap] Cloning SECQUOIA/QuIP into $repo_dir")
            run(`git clone --depth 1 https://github.com/SECQUOIA/QuIP.git $repo_dir`)
        end
        include(joinpath(repo_dir, "scripts", "notebook_bootstrap.jl"))
        return nothing
    end

    error("Could not locate scripts/notebook_bootstrap.jl from $(pwd()).")
end

load_notebook_bootstrap()

BOOTSTRAP = QuIPNotebookBootstrap.bootstrap_notebook("3-GAMA")

QUIP_REPO_DIR = BOOTSTRAP.repo_dir
JULIA_NOTEBOOKS_DIR = BOOTSTRAP.notebooks_dir
JULIA_PROJECT_DIR = BOOTSTRAP.project_dir
IN_COLAB = BOOTSTRAP.in_colab

SELECTED_GAMA_EXAMPLE = 4  # Change to 1, 2, 3, or 4, then rerun from this cell downward.

"""Return one of the class-reference GAMA examples."""
function load_gama_example(example_id::Integer)
    if example_id == 1
        return (
            id = 1,
            name = "Example 1: illustrative Graver augmentation",
            A = [1 1 1 1; 1 5 10 25],
            b = [21, 156],
            x0 = [1, 15, 3, 2],
            xl = zeros(Int, 4),
            xu = fill(15, 4),
            c = nothing,
            objective = :distance_to_5,
            graver_directions = [
                5 -6 0 1
                5 -9 4 0
                0 3 -4 1
                5 -3 -4 2
                5 0 -8 3
            ],
        )
    elseif example_id == 2
        return (
            id = 2,
            name = "Example 2: four-variable linear example",
            A = [1 1 1 1; 1 2 3 4],
            b = [10, 21],
            x0 = [1, 8, 0, 1],
            xl = zeros(Int, 4),
            xu = fill(21, 4),
            c = [0, 1, 0, 2],
            objective = :linear,
            graver_directions = [
                1 -2 1 0
                2 -3 0 1
                1 -1 -1 1
                0 1 -2 1
                1 0 -3 2
            ],
        )
    elseif example_id == 3
        return (
            id = 3,
            name = "Example 3: alternate four-variable linear example",
            A = [1 1 1 1; 0 1 2 3],
            b = [10, 15],
            x0 = [3, 0, 6, 1],
            xl = zeros(Int, 4),
            xu = fill(10, 4),
            c = [1, 3, 14, 17],
            objective = :linear,
            graver_directions = [
                2 -3 0 1
                1 -2 1 0
                1 -1 -1 1
                0 1 -2 1
                1 0 -3 2
            ],
        )
    elseif example_id == 4
        A = [
            1 1 1 1 1 1 1 1 1 1 0 1 0 1 0 1 0 1 1 1 0 1 0 1 0
            1 1 1 1 0 1 0 1 0 0 1 0 0 0 1 0 0 1 0 1 1 1 1 1 1
            0 1 0 0 0 1 0 1 0 1 1 0 1 1 0 1 1 0 0 1 0 0 1 1 1
            0 0 0 0 0 0 0 1 0 1 1 1 0 1 1 1 1 0 0 1 0 0 0 0 0
            0 1 1 1 1 1 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0 1 0 1 0
        ]
        n = size(A, 2)
        return (
            id = 4,
            name = "Example 4: saved portfolio instance",
            A = A,
            b = vec(ceil.(Int, sum(A; dims = 2) / 2)),
            x0 = [1, 1, 1, 1, -1, 1, 1, 1, 1, 1, 1, 1, -2, 1, 0, -1, 0, 1, -1, 1, -2, -2, 1, 1, 1],
            xl = fill(-2, n),
            xu = fill(2, n),
            c = nothing,
            objective = :portfolio,
            graver_directions = nothing,
        )
    else
        error("Unknown GAMA example $(example_id); choose 1, 2, 3, or 4.")
    end
end

gama_example = load_gama_example(SELECTED_GAMA_EXAMPLE)
A = gama_example.A
b = gama_example.b
x0 = gama_example.x0
xl = gama_example.xl
xu = gama_example.xu
c = gama_example.c
m, n = size(A)

println("Running $(gama_example.name) with A size $(size(A)) and $(length(x0)) variables.")

println("Initial point is feasible: ", A * x0 == b)
println("Variable bounds: [", minimum(xl), ", ", maximum(xu), "]")

if !@isdefined(DelimitedFiles)
    using DelimitedFiles
end

if !@isdefined(Random)
    using Random
end

if !@isdefined(LinearAlgebra)
    using LinearAlgebra
end

"""Return the path to a shared Notebook 3 data file."""
shared_data_path(filename::AbstractString) = joinpath(QUIP_REPO_DIR, "notebooks_data", filename)

"""
Return the fixed `μ` and `σ` vectors shared by both notebooks.

Returns
-------
Tuple{Vector{Float64}, Vector{Float64}}
    The expected-return vector `μ` and the volatility vector `σ`.
"""
function load_shared_coefficients(num_variables::Integer)
    coeff_path = shared_data_path("3-GAMA_example4_coefficients.csv")
    if !isfile(coeff_path)
        rng = MersenneTwister(1729)
        μ = rand(rng, num_variables)
        σ = rand(rng, num_variables) .* μ
        println("Shared coefficient file not found; generating deterministic coefficients with seed 1729.")
        return (μ, σ)
    end

    coeffs = readdlm(coeff_path, ',', Float64)
    size(coeffs) == (2, num_variables) || error("Expected coefficient array with shape (2, $(num_variables)), found $(size(coeffs)).")
    return (vec(coeffs[1, :]), vec(coeffs[2, :]))
end

ε = 0.01
if SELECTED_GAMA_EXAMPLE == 4
    μ, σ = load_shared_coefficients(n)
else
    μ = zeros(Float64, n)
    σ = zeros(Float64, n)
end

function f(x)
    if SELECTED_GAMA_EXAMPLE == 1
        return sum(abs.(x .- 5))
    elseif SELECTED_GAMA_EXAMPLE == 4
        return -dot(μ, x) + sqrt(((1 - ε) / ε) * dot(σ .^ 2, x .^ 2))
    else
        return dot(c, x)
    end
end

using DelimitedFiles
import lib4ti2_jll
import BinaryWrappers
using NPZ

const lib4ti2_bin = BinaryWrappers.@generate_wrappers(lib4ti2_jll)

"""Return `true` when the bundled 4ti2 `graver` executable is available."""
function has_graver()::Bool
    try
        return success(`$(lib4ti2_bin)/graver --help`)
    catch
        return false
    end
end

"""Print the help text for the bundled 4ti2 `graver` executable."""
function graver()
    run(`$(lib4ti2_bin)/graver --help`)

    return nothing
end

"""Run the bundled 4ti2 `graver` executable on the given project path."""
function graver(proj_path::AbstractString; silent::Bool = true)
    if silent
        run(`$(lib4ti2_bin)/graver -q $(proj_path)`)
    else
        run(`$(lib4ti2_bin)/graver $(proj_path)`)
    end

    return nothing
end

"""Write an integer matrix in the text format expected by 4ti2."""
function write_mat(path::AbstractString, A)
    m, n = size(A)

    open(path, "w") do io
        println(io, "$m $n")

        join(io, (join(@view(A[i, :]), " ") for i = 1:m), "\n")
    end

    return nothing
end

"""Read a matrix written in 4ti2's plain-text format."""
function read_mat(path::AbstractString, type = Int)
    m, n = parse.(Int, split(readline(path)))

    A = Matrix{type}(undef, m, n)

    open(path, "r") do io
        readline(io)

        for (i, line) in enumerate(eachline(io))
            A[i, :] .= parse.(Int, split(line))
        end
    end

    return A
end

"""
Return the Graver basis of `A` by calling the bundled 4ti2 `graver` executable.

Only the constraint matrix is written because the Graver basis depends on `A` alone;
the box bounds `xl` and `xu` are enforced later during augmentation.
"""
function compute_graver_basis_local(A)
    G = nothing

    mktempdir() do path
        proj_path = joinpath(path, "proj")

        write_mat("$(proj_path).mat", A)
        graver(proj_path)

        G = read_mat("$(proj_path).gra")
    end

    return G
end

"""Load the bundled `graver.npy` file when the 4ti2 executable is unavailable."""
function download_graver_basis()
    npy_path = joinpath(JULIA_NOTEBOOKS_DIR, "graver.npy")
    G = NPZ.npyread(npy_path)
    return Array{Int}(G)
end

"""Return the Graver basis of `A`, preferring the bundled 4ti2 binary when available."""
function graver_basis(A)
    if has_graver()
        return compute_graver_basis_local(A)
    else
        return download_graver_basis()
    end
end

function graver_basis_for_selected_example(A)
    if gama_example.graver_directions !== nothing
        println("Using $(size(gama_example.graver_directions, 1)) built-in Graver directions for Example $(SELECTED_GAMA_EXAMPLE).")
        return gama_example.graver_directions
    end

    return graver_basis(A)
end

G = graver_basis_for_selected_example(A)

# Define rules to choose augmentation element, either the best one (argmin) or the first one that is found
function argmin_rule(a)
    local i, y, α = (nothing, NaN, Inf)

    for (_i, (_y, _α)) in enumerate(a)
        if _α < α
            i, y, α = (_i, _y, _α)
        end
    end

    return (i, (y, α))
end

function greedy_rule(a)
    local i, y, α = (nothing, NaN, Inf)

    for (_i, (_y, _α)) in enumerate(a)
        if _α ≉ 0.0
            return (_i, (_y, _α))
        else
            i, y, α = (_i, _y, _α)
        end
    end

    return (i, (y, α))
end

# Bisection rules for finding best step size
"""Compute the best integer step size along a Graver direction within the box bounds."""
function bisection_rule(f::Function, g, x, xl = nothing, xu = nothing, laststep = nothing)
    if g == laststep
        return (f(x), 0)
    end

    if xl === nothing
        xl = zeros(length(x))
    end

    if xu === nothing
        xu = fill(2 * maximum(x), length(x))
    end

    u = maximum(xu) - minimum(xl)
    l = -u

    for (i, gi) in enumerate(g)
        if gi >= 1
            u = min(u, floor(Int, (xu[i] - x[i]) / gi))
            l = max(l, ceil(Int, (xl[i] - x[i]) / gi))
        elseif gi <= -1
            u = min(u, floor(Int, (xl[i] - x[i]) / gi))
            l = max(l, ceil(Int, (xu[i] - x[i]) / gi))
        end
    end

    α = u

    while u - l > 1
        α = ifelse(f(x + l * g) < f(x + u * g), l, u)

        p1 = floor(Int, (l + u) / 2) - 1
        p2 = floor(Int, (l + u) / 2)
        p3 = floor(Int, (l + u) / 2) + 1

        if f(x + p1 * g) < f(x + p2 * g)
            u = floor(Int, (l + u) / 2)
        elseif f(x + p3 * g) < f(x + p2 * g)
            l = floor(Int, (l + u) / 2) + 1
        else
            α = p2
            break
        end
    end

    if f(x + l * g) < f(x + u * g) && f(x + l * g) < f(x + α * g)
        α = l
    elseif f(x + u * g) < f(x + α * g)
        α = u
    end

    return (f(x + α * g), α)
end

# We can just have a single step move (works well with greedy approach)
"""Evaluate whether a single forward or backward step improves the objective."""
function single_move_rule(f::Function, g, x, xl = nothing, xu = nothing, laststep = nothing)
    if xl === nothing
        xl = zeros(length(x))
    end

    if xu === nothing
        xu = fill(2 * maximum(x), length(x))
    end

    α = 0

    if all(x + g .<= xu) && all(x + g .>= xl)
        if f(x + g) < f(x)
            α = 1
        end
    elseif all(x - g .<= xu) && all(x - g .>= xl)
        if f(x - g) < f(x) && f(x - g) < f(x + g)
            α = -1
        end
    end

    return (f(x + α * g), α)
end

"""Compute the Graver basis for `A` and run the augmentation experiment from `x0`."""
function augmentation(
    f::Function,
    A,
    b,
    x0,
    xl,
    xu;
    silent::Bool = true,
    max_iter::Integer = 1_000,
    step_rule::Function = bisection_rule,
    choice_rule::Function = argmin_rule,
)
    G = graver_basis_for_selected_example(A)

    return augmentation(
        f,
        G,
        A,
        b,
        x0,
        xl,
        xu;
        silent,
        max_iter,
        step_rule,
        choice_rule,
    )
end

"""Apply the selected augmentation rule to a fixed Graver basis until convergence."""
function augmentation(
    f::Function,
    G,
    A,
    b,
    x0,
    xl,
    xu;
    silent::Bool = true,
    max_iter::Integer = 1_000,
    step_rule::Function = bisection_rule,
    choice_rule::Function = argmin_rule,
)
    # Constraints definition
    isfeasible = (x) -> (A * x == b)

    Δ = 1
    n = length(x0)
    k = 1

    if !silent
        println("Initial point: $(x0)")
        println("Objective function: $(f(x0))")
    end

    y = nothing
    x = copy(x0)
    δ = Vector{Int}(undef, n)
    gprev = Vector{Int}(undef, n)

    while Δ != 0 && k < max_iter
        i, (y, Δ) = choice_rule((step_rule(f, g, x, xl, xu, gprev) for g in eachrow(G)))

        δ .= @view(G[i, :]) .* Δ
        x .+= δ
        gprev .= @view(G[i, :])

        if !silent
            println("Iteration ", k)
            println(i, (y, Δ))
            println("Augmentation direction:", gprev)
            println("Distanced moved:", Δ)
            println("Step taken:", δ)
            println("Objective function:", y)
            println(f(x))
            println("Current point:", x)
            println("Are constraints satisfied?", isfeasible(x))
        end

        k += 1
    end

    return (k, y, x)
end

println("Best-augmentation: Choosing among the best step that each element of G can do (via bisection), the one that reduces the most the objective")

iter, f_obj, xf = @time augmentation(
    f, G, A, b, x0, xl, xu;
    step_rule = bisection_rule, choice_rule = argmin_rule
)

println("$(iter), iterations")
println("solution: $(xf)")
println("objective: $(f_obj)")

println("Greedy-best-augmentation: Choosing among the best step that each element of G can do (via bisection), the first one encountered that reduces the objective")

iter,f_obj,xf = @time augmentation(
    f, G, A, b, x0, xl, xu;
    step_rule = bisection_rule, choice_rule = greedy_rule,
)

println("$(iter), iterations")
println("solution: $(xf)")
println("objective: $(f_obj)")

println("Greedy-augmentation: Choosing among the first element of G that with a single step reduces the objective")

iter,f_obj,xf = @time augmentation(
    f, G, A, b, x0, xl, xu;
    step_rule = single_move_rule, choice_rule = greedy_rule,
)

println("$(iter), iterations")
println("solution: $(xf)")
println("objective: $(f_obj)")

if !@isdefined(JuMP)
    using JuMP
end

if !@isdefined(DWave)
    using DWave
end

if !@isdefined(LinearAlgebra)
    using LinearAlgebra
end

"""Return the committed Example 4 feasible starts when available."""
function load_precomputed_feasible_starts()
    SELECTED_GAMA_EXAMPLE == 4 || return nothing

    feasible_path = shared_data_path("3-GAMA_example4_feasible_starts.csv")
    if !isfile(feasible_path)
        return nothing
    end

    X = readdlm(feasible_path, ',', Int)
    println("Loaded $(size(X, 1)) precomputed feasible solutions from 3-GAMA_example4_feasible_starts.csv.")
    return [vec(X[i, :]) for i in axes(X, 1)]
end

"""Return Example 4 starts, or use the documented start for smaller examples."""
function get_default_feasible_starts(A, b; num_reads = 20)
    if SELECTED_GAMA_EXAMPLE != 4
        println("Using the documented x0 as the feasible start for Example $(SELECTED_GAMA_EXAMPLE).")
        return [copy(x0)]
    end

    X = load_precomputed_feasible_starts()
    if X !== nothing
        return X
    end

    println("Shared feasible-start file not found; recomputing feasible solutions with simulated annealing.")
    return get_feasible(A, b; num_reads = num_reads)
end

"""Return the saved Example 4 Graver order, or the natural order for small examples."""
function load_graver_order(num_directions::Int)
    if SELECTED_GAMA_EXAMPLE != 4
        return collect(1:num_directions)
    end

    order_path = shared_data_path("3-GAMA_example4_graver_order.csv")
    if isfile(order_path)
        order = vec(readdlm(order_path, ',', Int)) .+ 1
        length(order) == num_directions || error("Expected $(num_directions) Graver indices, found $(length(order)).")
        return order
    end

    rng = MersenneTwister(271828)
    println("Shared Graver-order file not found; generating a deterministic permutation with seed 271828.")
    return randperm(rng, num_directions)
end

"""
Sample binary feasible starting points by minimizing the QUBO residual.

Arguments
---------
A, b:
    Constraint matrix and right-hand side.
num_reads:
    Number of reads requested from the annealer.

Returns
-------
Vector{Vector{Int}}
    Binary feasible starting points with one vector per sample.
"""
function get_feasible(A, b; num_reads = 1_000)
    m, n = size(A)

    Q = A' * A - 2 * diagm(A' * b)
    β = b' * b

    # Define Binary Quadratic Model using new DWave.jl
    model = Model(DWave.Neal.Optimizer)

    @variable(model, x[1:n], Bin)

    @objective(model, Min, x' * Q * x + β)

    set_optimizer_attribute(model, "num_reads", num_reads)

    optimize!(model)

    X = [
        round.(Int, value.(x; result = i))
        for i = 1:result_count(model)
        if objective_value(model; result = i) ≈ 0.0
    ]

    println("$(length(X)) feasible solutions found.")

    return X
end

X_feas = get_default_feasible_starts(A, b; num_reads = 20)
G_order = load_graver_order(size(G, 1));

"""
Run the full-basis augmentation experiment for each feasible starting point in `X`.

Arguments
---------
f, G, A, b, X, xl, xu:
    Objective function, Graver basis, linear constraints, feasible starts, and box bounds.

Returns
-------
Tuple{Vector{Float64}, Vector{Int}, Vector{Float64}}
    Final objectives, iteration counts, and runtimes for each feasible start.
"""
function augmentation_experiment(f, G, A, b, X, xl, xu)
    K = length(X)

    Y = Vector{Float64}(undef, K)
    I = Vector{Int}(undef, K)
    T = Vector{Float64}(undef, K)

    for (i, x) in enumerate(X)
        result = @timed augmentation(
            f, G, A, b, x, xl, xu;
            step_rule   = single_move_rule,
            choice_rule = greedy_rule,
        )

        num_iter, obj_val, _ = result.value

        Y[i] = obj_val
        I[i] = num_iter
        T[i] = result.time
    end

    return (Y, I, T)
end

Y_feas = f.(X_feas)

Y_aug, I_aug, T_aug = augmentation_experiment(f, G, A, b, X_feas, xl, xu);

if !@isdefined(Measures)
    using Measures
end

if !@isdefined(Random)
    using Random
end

if !@isdefined(Plots)
    using Plots
end

if !@isdefined(StatsBase)
    using StatsBase
end

if !@isdefined(StatsPlots)
    using StatsPlots
end

if !@isdefined(Printf)
    using Printf
end

"""
Plot the initial and augmented objectives together with the iteration counts.

Arguments
---------
Y_feas, Y_aug, I_aug:
    Initial objectives, augmented objectives, and iteration counts.
experiment_name:
    Title shown on the plot.

Returns
-------
Plots.Plot
    The combined plot object.
"""
function plot_augmentation(Y_feas, Y_aug, I_aug; experiment_name = "Augmentation")
    plt = plot(;
        plot_title = experiment_name,
        top_margin = 6mm,
        bottom_margin = 15mm,
    )

    plot!(
        plt, Y_feas;
        label       = "Initial",
        color       = :blue,
        markershape = :square,
    )
    plot!(
        plt, Y_aug;
        label       = "Augmented",
        color       = :orange,
        markershape = :square,
        ylabel      = "Objective Value",
        xlabel      = "Solutions",
        legend      = (0.1, -0.25),
    )

    scatter!(
        twinx(plt), I_aug;
        label       = "Iterations",
        color       = :green,
        legend      = (0.85, -0.25),
        markershape = :diamond,
        ylabel      = "Iterations",
        ylims       = (0, maximum(I_aug) + 1),
    )

    return plt
end

plot_augmentation(Y_feas, Y_aug, I_aug; experiment_name = "Complete-basis greedy ($(size(G, 1)) directions)")

"""
Run the augmentation experiment on a shared subset of the Graver basis.

Arguments
---------
f, G, A, b, X, xl, xu:
    Objective function, Graver basis, linear constraints, feasible starts, and box bounds.
order:
    Shared ordering of the Graver directions used to build nested subsets.
num_samples:
    Number of Graver directions kept in the subset.

Returns
-------
Tuple{Vector{Float64}, Vector{Int}, Vector{Float64}}
    Final objectives, iteration counts, and runtimes for each feasible start.
"""
function partial_augmentation_experiment(f, G, A, b, X, xl, xu, order; num_samples = 10)
    num_samples = min(num_samples, length(order), size(G, 1))
    K = length(X)

    Y = Vector{Float64}(undef, K)
    I = Vector{Int}(undef, K)
    T = Vector{Float64}(undef, K)

    index = order[1:num_samples]
    G_idx = G[index, :]

    for (i, x) in enumerate(X)
        result = @timed augmentation(
            f, G_idx, A, b, x, xl, xu;
            step_rule   = single_move_rule,
            choice_rule = greedy_rule,
        )

        num_iter, obj_val, _ = result.value

        Y[i] = obj_val
        I[i] = num_iter
        T[i] = result.time
    end

    return (Y, I, T)
end

num_partial_directions = min(10, size(G, 1))
Y_feas = f.(X_feas)

Y_paug, I_paug, T_paug = partial_augmentation_experiment(f, G, A, b, X_feas, xl, xu, G_order; num_samples = num_partial_directions);

plot_augmentation(Y_feas, Y_paug, I_paug; experiment_name = "Partial-basis greedy ($(num_partial_directions) directions)")

function log_ticks_for(values; include_zero_floor = nothing)
    positive = [Float64(v) for v in vec(values) if Float64(v) > 0]
    if isempty(positive)
        ticks = [1.0]
    else
        tick_lo = floor(Int, log10(minimum(positive)))
        tick_hi = ceil(Int, log10(maximum(positive)))
        ticks = 10.0 .^ collect(tick_lo:tick_hi)
    end

    if include_zero_floor !== nothing
        ticks = sort(unique(vcat([Float64(include_zero_floor)], ticks)))
    end

    labels = [include_zero_floor !== nothing && isapprox(t, Float64(include_zero_floor); atol = max(Float64(include_zero_floor) / 10, eps(Float64)), rtol = 0.0) ? "0" : @sprintf("%.1e", t) for t in ticks]
    return ticks, labels
end

"""
Plot the runtime comparison between the full and partial augmentation experiments.

Arguments
---------
T_aug, T_paug:
    Full-basis and partial-basis runtimes.
partial_label:
    Legend label used for the partial-basis series.

Returns
-------
Plots.Plot
    The runtime comparison plot.
"""
function plot_augmentation_runtime(T_aug, T_paug; partial_label = "10 sampled Graver directions")
    ticks, labels = log_ticks_for(vcat(T_aug, T_paug))
    plt = plot(;
        plot_title = "Augmentation Runtime",
        xlabel     = "Solutions",
        ylabel     = "Runtime (s)",
        legend     = :outertopright,
        right_margin = 18mm,
        size       = (760, 400),
        yscale     = :log10,
        yticks     = (ticks, labels),
    )

    plot!(T_aug; label = "Complete basis", markershape = :diamond)
    plot!(T_paug; label = partial_label, markershape = :diamond)

    return plt
end

plot_augmentation_runtime(T_aug, T_paug; partial_label = "$(num_partial_directions) sampled Graver directions")

"""
Run the partial-basis experiment over several fractions of the full Graver basis.

Arguments
---------
f, G, A, b, X, xl, xu:
    Objective function, Graver basis, linear constraints, feasible starts, and box bounds.
order:
    Shared ordering of the Graver directions used to build nested subsets.

Returns
-------
Tuple{Matrix{Float64}, Matrix{Int}, Matrix{Float64}}
    Final objectives, iteration counts, and runtimes for each sample-size fraction.
"""
function multiple_partial_augmentation_experiment(f, G, A, b, X, xl, xu, order)
    N = 10 # Discretization of the fractions of Graver considered
    M = length(X)
    K = size(G, 1)

    Y = Matrix{Float64}(undef, N, M)
    I = Matrix{Int}(undef, N, M)
    T = Matrix{Float64}(undef, N, M)

    for j in 1:N
        num_samples = max(1, trunc(Int, K / N * j))
        Y_p, I_p, T_p = partial_augmentation_experiment(
            f, G, A, b, X, xl, xu, order; num_samples = num_samples
        )

        Y[j, :] .= Y_p
        I[j, :] .= I_p
        T[j, :] .= T_p
    end

    return (Y, I, T)
end

Y_mpaug, I_mpaug, T_mpaug = multiple_partial_augmentation_experiment(f, G, A, b, X_feas, xl, xu, G_order);

"""Return log-safe objective gaps together with the floor used to display zero-gap entries."""
function lift_zero_gaps(Y, global_minimum)
    Y_gap = Y .- global_minimum
    positive = Y_gap[Y_gap .> 0]
    zero_floor = isempty(positive) ? eps(Float64) : minimum(positive) / 10
    adjusted = map(y -> y > 0 ? y : zero_floor, Y_gap)
    return adjusted, zero_floor
end

"""
Plot the log-scaled objective gap to the best full-basis result across sample sizes.

Arguments
---------
Y_feas, Y_mpaug, global_minimum:
    Initial objectives, partial-basis objectives, and the best full-basis objective.

Returns
-------
Plots.Plot
    The boxplot object with zero-gap entries lifted for the log axis.
"""
function plot_multiple_partial_augmentation(Y_feas, Y_mpaug, global_minimum)
    Y = hcat(Y_feas, Y_mpaug')
    Y_gap, zero_floor = lift_zero_gaps(Y, global_minimum)
    X = permutedims(["Initial"; ["$(10i)% |G|" for i = 1:size(Y_mpaug, 1)]])

    ticks, labels = log_ticks_for(Y_gap; include_zero_floor = zero_floor)

    plt = boxplot(X, Y_gap;
        plot_title = "Multiple partial augmentation samples",
        size       = (800, 400),
        margin     = 5mm,
        ylabel     = "Objective gap to best full-basis result",
        xlabel     = "Sample Sizes",
        yscale     = :log10,
        yticks     = (ticks, labels),
        legend     = false,
    )

    return plt
end

plot_multiple_partial_augmentation(Y_feas, Y_mpaug, minimum(Y_aug))