Replay

src/BB/log.jl

@@ -31,9 +31,9 @@ function print_progress(
     node::Node;
     time_elapsed::Float64,
     print_interval::Int,
-    primal_update::Bool,
+    primal_update::Bool
 )::Nothing
-    if (pool.processed % print_interval == 0) || isempty(pool.pending) || primal_update
+    if (pool.processed % print_interval == 0) || isempty(pool.pending)
         if isempty(node.branch_vars)
             branch_var_name = "---"
             branch_lb = "---"

src/BB/nodepool.jl

@@ -8,10 +8,10 @@ import Base.Threads: threadid
 
 function take(
     pool::NodePool;
-    suggestions::Array{Node} = [],
+    suggestions::Array{Node}=[],
     time_remaining::Float64,
     gap_limit::Float64,
-    node_limit::Int,
+    node_limit::Int
 )::Union{Symbol,Node}
     t = threadid()
     lock(pool.lock) do
@@ -53,8 +53,8 @@ function offer(
     pool::NodePool;
     parent_node::Union{Nothing,Node},
     child_nodes::Vector{Node},
-    time_elapsed::Float64 = 0.0,
-    print_interval::Int = 100,
+    time_elapsed::Float64=0.0,
+    print_interval::Int=100
 )::Nothing
     lock(pool.lock) do
         primal_update = false
@@ -101,15 +101,16 @@ function offer(
             # Update branching variable history
             branch_var = child_nodes[1].branch_vars[end]
             offset = findfirst(isequal(branch_var), parent_node.fractional_variables)
+            if offset !== nothing
                 x = parent_node.fractional_values[offset]
                 obj_change_up = child_nodes[1].obj - parent_node.obj
                 obj_change_down = child_nodes[2].obj - parent_node.obj
                 _update_var_history(
-                pool = pool,
-                var = branch_var,
-                x = x,
-                obj_change_down = obj_change_down,
-                obj_change_up = obj_change_up,
+                    pool=pool,
+                    var=branch_var,
+                    x=x,
+                    obj_change_down=obj_change_down,
+                    obj_change_up=obj_change_up,
                 )
                 # Update global history
                 pool.history.avg_pseudocost_up =
@@ -117,14 +118,15 @@ function offer(
                 pool.history.avg_pseudocost_down =
                     mean(vh.pseudocost_down for vh in values(pool.var_history))
             end
+        end
 
         for node in child_nodes
             print_progress(
                 pool,
                 node,
-                time_elapsed = time_elapsed,
-                print_interval = print_interval,
-                primal_update = isfinite(node.obj) && isempty(node.fractional_variables),
+                time_elapsed=time_elapsed,
+                print_interval=print_interval,
+                primal_update=isfinite(node.obj) && isempty(node.fractional_variables),
             )
         end
     end
@@ -136,7 +138,7 @@ function _update_var_history(;
     var::Variable,
     x::Float64,
     obj_change_down::Float64,
-    obj_change_up::Float64,
+    obj_change_up::Float64
 )::Nothing
     # Create new history entry
     if var ∉ keys(pool.var_history)

src/BB/optimize.jl

@@ -10,16 +10,22 @@ import ..H5File
 
 function solve!(
     mip::MIP;
-    time_limit::Float64 = Inf,
-    node_limit::Int = typemax(Int),
-    gap_limit::Float64 = 1e-4,
-    print_interval::Int = 5,
-    initial_primal_bound::Float64 = Inf,
-    branch_rule::VariableBranchingRule = ReliabilityBranching(),
-    enable_plunging = true,
-)::NodePool
+    time_limit::Float64=Inf,
+    node_limit::Int=typemax(Int),
+    gap_limit::Float64=1e-4,
+    print_interval::Int=5,
+    initial_primal_bound::Float64=Inf,
+    branch_rule::VariableBranchingRule=ReliabilityBranching(),
+    enable_plunging=true,
+    replay=nothing
+)::Tuple{NodePool,ReplayInfo}
+
+    if replay === nothing
+        replay = ReplayInfo()
+    end
+
     time_initial = time()
-    pool = NodePool(mip = mip)
+    pool = NodePool(mip=mip, next_index=replay.next_index)
     pool.primal_bound = initial_primal_bound
 
     root_node = _create_node(mip)
@@ -34,9 +40,9 @@ function solve!(
 
     offer(
         pool,
-        parent_node = nothing,
-        child_nodes = [root_node],
-        print_interval = print_interval,
+        parent_node=nothing,
+        child_nodes=[root_node],
+        print_interval=print_interval,
     )
     @threads for t = 1:nthreads()
         child_one, child_zero, suggestions = nothing, nothing, Node[]
@@ -47,10 +53,10 @@ function solve!(
             end
             node = take(
                 pool,
-                suggestions = suggestions,
-                time_remaining = time_limit - time_elapsed,
-                node_limit = node_limit,
-                gap_limit = gap_limit,
+                suggestions=suggestions,
+                time_remaining=time_limit - time_elapsed,
+                node_limit=node_limit,
+                gap_limit=gap_limit,
             )
             if node == :END
                 break
@@ -64,10 +70,25 @@ function solve!(
                 @assert status == :Optimal
                 _unset_node_bounds(node)
 
+                if node.index in keys(replay.node_decisions)
+                    decision = replay.node_decisions[node.index]
+                    ids = decision.ids
+                    branch_var = decision.branch_var
+                    var_value = decision.var_value
+                else
                     # Find branching variable
                     ids = generate_indices(pool, 2)
                     branch_var = find_branching_var(branch_rule, node, pool)
 
+                    # Query current fractional value
+                    offset = findfirst(isequal(branch_var), node.fractional_variables)
+                    var_value = node.fractional_values[offset]
+
+                    # Update replay
+                    decision = ReplayNodeDecision(; branch_var, var_value, ids)
+                    replay.node_decisions[node.index] = decision
+                end
+
                 # Find current variable lower and upper bounds
                 offset = findfirst(isequal(branch_var), mip.int_vars)
                 var_lb = mip.int_vars_lb[offset]
@@ -79,46 +100,43 @@ function solve!(
                     end
                 end
 
-                # Query current fractional value
-                offset = findfirst(isequal(branch_var), node.fractional_variables)
-                var_value = node.fractional_values[offset]
-
                 child_zero = _create_node(
                     mip,
-                    index = ids[2],
-                    parent = node,
-                    branch_var = branch_var,
-                    branch_var_lb = var_lb,
-                    branch_var_ub = floor(var_value),
+                    index=ids[2],
+                    parent=node,
+                    branch_var=branch_var,
+                    branch_var_lb=var_lb,
+                    branch_var_ub=floor(var_value),
                 )
                 child_one = _create_node(
                     mip,
-                    index = ids[1],
-                    parent = node,
-                    branch_var = branch_var,
-                    branch_var_lb = ceil(var_value),
-                    branch_var_ub = var_ub,
+                    index=ids[1],
+                    parent=node,
+                    branch_var=branch_var,
+                    branch_var_lb=ceil(var_value),
+                    branch_var_ub=var_ub,
                 )
                 offer(
                     pool,
-                    parent_node = node,
-                    child_nodes = [child_one, child_zero],
-                    time_elapsed = time_elapsed,
-                    print_interval = print_interval,
+                    parent_node=node,
+                    child_nodes=[child_one, child_zero],
+                    time_elapsed=time_elapsed,
+                    print_interval=print_interval,
                 )
             end
         end
     end
-    return pool
+    replay.next_index = pool.next_index
+    return pool, replay
 end
 
 function _create_node(
     mip;
-    index::Int = 0,
-    parent::Union{Nothing,Node} = nothing,
-    branch_var::Union{Nothing,Variable} = nothing,
-    branch_var_lb::Union{Nothing,Float64} = nothing,
-    branch_var_ub::Union{Nothing,Float64} = nothing,
+    index::Int=0,
+    parent::Union{Nothing,Node}=nothing,
+    branch_var::Union{Nothing,Variable}=nothing,
+    branch_var_lb::Union{Nothing,Float64}=nothing,
+    branch_var_ub::Union{Nothing,Float64}=nothing
 )::Node
     if parent === nothing
         branch_vars = Variable[]

src/BB/structs.jl

@@ -72,3 +72,14 @@ Base.@kwdef mutable struct NodePool
     history::History = History()
     var_history::Dict{Variable,VariableHistory} = Dict()
 end
+
+Base.@kwdef struct ReplayNodeDecision
+    branch_var
+    var_value
+    ids
+end
+
+Base.@kwdef mutable struct ReplayInfo
+    node_decisions::Dict{Int,ReplayNodeDecision} = Dict()
+    next_index::Int = 1
+end

src/solvers/jump.jl

@@ -36,13 +36,14 @@ function _add_constrs(
 end
 
 function _extract_after_load(model::JuMP.Model, h5)
+    @info "_extract_after_load"
     if JuMP.objective_sense(model) == MOI.MIN_SENSE
         h5.put_scalar("static_sense", "min")
     else
         h5.put_scalar("static_sense", "max")
     end
-    _extract_after_load_vars(model, h5)
-    _extract_after_load_constrs(model, h5)
+    @time _extract_after_load_vars(model, h5)
+    @time _extract_after_load_constrs(model, h5)
 end
 
 function _extract_after_load_vars(model::JuMP.Model, h5)
@@ -117,10 +118,11 @@ function _extract_after_load_constrs(model::JuMP.Model, h5)
 end
 
 function _extract_after_lp(model::JuMP.Model, h5)
+    @info "_extract_after_lp"
     h5.put_scalar("lp_wallclock_time", solve_time(model))
     h5.put_scalar("lp_obj_value", objective_value(model))
-    _extract_after_lp_vars(model, h5)
-    _extract_after_lp_constrs(model, h5)
+    @time _extract_after_lp_vars(model, h5)
+    @time _extract_after_lp_constrs(model, h5)
 end
 
 function _extract_after_lp_vars(model::JuMP.Model, h5)
@@ -146,46 +148,46 @@ function _extract_after_lp_vars(model::JuMP.Model, h5)
     end
     h5.put_array("lp_var_basis_status", to_str_array(basis_status))
 
-    # Sensitivity analysis
-    obj_coeffs = h5.get_array("static_var_obj_coeffs")
-    sensitivity_report = lp_sensitivity_report(model)
-    sa_obj_down, sa_obj_up = Float64[], Float64[]
-    sa_lb_down, sa_lb_up = Float64[], Float64[]
-    sa_ub_down, sa_ub_up = Float64[], Float64[]
-    for (i, v) in enumerate(vars)
-        # Objective function
-        (delta_down, delta_up) = sensitivity_report[v]
-        push!(sa_obj_down, delta_down + obj_coeffs[i])
-        push!(sa_obj_up, delta_up + obj_coeffs[i])
-
-        # Lower bound
-        if has_lower_bound(v)
-            constr = LowerBoundRef(v)
-            (delta_down, delta_up) = sensitivity_report[constr]
-            push!(sa_lb_down, lower_bound(v) + delta_down)
-            push!(sa_lb_up, lower_bound(v) + delta_up)
-        else
-            push!(sa_lb_down, -Inf)
-            push!(sa_lb_up, -Inf)
-        end
-
-        # Upper bound
-        if has_upper_bound(v)
-            constr = JuMP.UpperBoundRef(v)
-            (delta_down, delta_up) = sensitivity_report[constr]
-            push!(sa_ub_down, upper_bound(v) + delta_down)
-            push!(sa_ub_up, upper_bound(v) + delta_up)
-        else
-            push!(sa_ub_down, Inf)
-            push!(sa_ub_up, Inf)
-        end
-    end
-    h5.put_array("lp_var_sa_obj_up", sa_obj_up)
-    h5.put_array("lp_var_sa_obj_down", sa_obj_down)
-    h5.put_array("lp_var_sa_ub_up", sa_ub_up)
-    h5.put_array("lp_var_sa_ub_down", sa_ub_down)
-    h5.put_array("lp_var_sa_lb_up", sa_lb_up)
-    h5.put_array("lp_var_sa_lb_down", sa_lb_down)
+    # # Sensitivity analysis
+    # obj_coeffs = h5.get_array("static_var_obj_coeffs")
+    # sensitivity_report = lp_sensitivity_report(model)
+    # sa_obj_down, sa_obj_up = Float64[], Float64[]
+    # sa_lb_down, sa_lb_up = Float64[], Float64[]
+    # sa_ub_down, sa_ub_up = Float64[], Float64[]
+    # for (i, v) in enumerate(vars)
+    #     # Objective function
+    #     (delta_down, delta_up) = sensitivity_report[v]
+    #     push!(sa_obj_down, delta_down + obj_coeffs[i])
+    #     push!(sa_obj_up, delta_up + obj_coeffs[i])
+
+    #     # Lower bound
+    #     if has_lower_bound(v)
+    #         constr = LowerBoundRef(v)
+    #         (delta_down, delta_up) = sensitivity_report[constr]
+    #         push!(sa_lb_down, lower_bound(v) + delta_down)
+    #         push!(sa_lb_up, lower_bound(v) + delta_up)
+    #     else
+    #         push!(sa_lb_down, -Inf)
+    #         push!(sa_lb_up, -Inf)
+    #     end
+
+    #     # Upper bound
+    #     if has_upper_bound(v)
+    #         constr = JuMP.UpperBoundRef(v)
+    #         (delta_down, delta_up) = sensitivity_report[constr]
+    #         push!(sa_ub_down, upper_bound(v) + delta_down)
+    #         push!(sa_ub_up, upper_bound(v) + delta_up)
+    #     else
+    #         push!(sa_ub_down, Inf)
+    #         push!(sa_ub_up, Inf)
+    #     end
+    # end
+    # h5.put_array("lp_var_sa_obj_up", sa_obj_up)
+    # h5.put_array("lp_var_sa_obj_down", sa_obj_down)
+    # h5.put_array("lp_var_sa_ub_up", sa_ub_up)
+    # h5.put_array("lp_var_sa_ub_down", sa_ub_down)
+    # h5.put_array("lp_var_sa_lb_up", sa_lb_up)
+    # h5.put_array("lp_var_sa_lb_down", sa_lb_down)
 end
 
 
@@ -201,7 +203,7 @@ function _extract_after_lp_constrs(model::JuMP.Model, h5)
     duals = Float64[]
     basis_status = []
     constr_idx = 1
-    sensitivity_report = lp_sensitivity_report(model)
+    # sensitivity_report = lp_sensitivity_report(model)
     for (ftype, stype) in JuMP.list_of_constraint_types(model)
         for constr in JuMP.all_constraints(model, ftype, stype)
             length(JuMP.name(constr)) > 0 || continue
@@ -219,21 +221,22 @@ function _extract_after_lp_constrs(model::JuMP.Model, h5)
                 error("Unknown basis status: $b")
             end
 
-            # Sensitivity analysis
-            (delta_down, delta_up) = sensitivity_report[constr]
-            push!(sa_rhs_down, rhs[constr_idx] + delta_down)
-            push!(sa_rhs_up, rhs[constr_idx] + delta_up)
+            # # Sensitivity analysis
+            # (delta_down, delta_up) = sensitivity_report[constr]
+            # push!(sa_rhs_down, rhs[constr_idx] + delta_down)
+            # push!(sa_rhs_up, rhs[constr_idx] + delta_up)
 
             constr_idx += 1
         end
     end
     h5.put_array("lp_constr_dual_values", duals)
     h5.put_array("lp_constr_basis_status", to_str_array(basis_status))
-    h5.put_array("lp_constr_sa_rhs_up", sa_rhs_up)
-    h5.put_array("lp_constr_sa_rhs_down", sa_rhs_down)
+    # h5.put_array("lp_constr_sa_rhs_up", sa_rhs_up)
+    # h5.put_array("lp_constr_sa_rhs_down", sa_rhs_down)
 end
 
 function _extract_after_mip(model::JuMP.Model, h5)
+    @info "_extract_after_mip"
     h5.put_scalar("mip_obj_value", objective_value(model))
     h5.put_scalar("mip_obj_bound", objective_bound(model))
     h5.put_scalar("mip_wallclock_time", solve_time(model))
@@ -254,11 +257,14 @@ end
 function _fix_variables(model::JuMP.Model, var_names, var_values, stats)
     vars = [variable_by_name(model, v) for v in var_names]
     for (i, var) in enumerate(vars)
-        fix(var, var_values[i], force = true)
+        if isfinite(var_values[i])
+            fix(var, var_values[i], force=true)
+        end
     end
 end
 
 function _optimize(model::JuMP.Model)
+    @info "_optimize"
     optimize!(model)
     flush(stdout)
     Libc.flush_cstdio()
@@ -269,7 +275,7 @@ function _relax(model::JuMP.Model)
     relax_integrality(relaxed)
     # FIXME: Remove hardcoded optimizer
     set_optimizer(relaxed, HiGHS.Optimizer)
-    set_silent(relaxed)
+    # set_silent(relaxed)
     return relaxed
 end
 
@@ -303,7 +309,7 @@ function __init_solvers_jump__()
             constrs_lhs,
             constrs_sense,
             constrs_rhs,
-            stats = nothing,
+            stats=nothing,
         ) = _add_constrs(
             self.inner,
             from_str_array(var_names),
@@ -319,14 +325,14 @@ function __init_solvers_jump__()
 
         extract_after_mip(self, h5) = _extract_after_mip(self.inner, h5)
 
-        fix_variables(self, var_names, var_values, stats = nothing) =
+        fix_variables(self, var_names, var_values, stats=nothing) =
             _fix_variables(self.inner, from_str_array(var_names), var_values, stats)
 
         optimize(self) = _optimize(self.inner)
 
         relax(self) = Class(_relax(self.inner))
 
-        set_warm_starts(self, var_names, var_values, stats = nothing) =
+        set_warm_starts(self, var_names, var_values, stats=nothing) =
             _set_warm_starts(self.inner, from_str_array(var_names), var_values, stats)
 
         write(self, filename) = _write(self.inner, filename)

test/Project.toml

@@ -4,7 +4,9 @@ authors = ["Alinson S. Xavier <git@axavier.org>"]
 version = "0.1.0"
 
 [deps]
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"

test/src/BB/test_bb.jl

@@ -10,9 +10,12 @@ using Test
 using MIPLearn.BB
 using MIPLearn
 
+using CSV
+using DataFrames
+
 basepath = @__DIR__
 
-function bb_run(optimizer_name, optimizer; large = true)
+function bb_run(optimizer_name, optimizer; large=true)
     @testset "Solve ($optimizer_name)" begin
         @testset "interface" begin
             filename = "$FIXTURES/danoint.mps.gz"
@@ -25,7 +28,7 @@ function bb_run(optimizer_name, optimizer; large = true)
 
             status, obj = BB.solve_relaxation!(mip)
             @test status == :Optimal
-            @test round(obj, digits = 6) == 62.637280
+            @test round(obj, digits=6) == 62.637280
 
             @test BB.name(mip, mip.int_vars[1]) == "xab"
             @test BB.name(mip, mip.int_vars[2]) == "xac"
@@ -35,26 +38,26 @@ function bb_run(optimizer_name, optimizer; large = true)
             @test mip.int_vars_ub[1] == 1.0
 
             vals = BB.values(mip, mip.int_vars)
-            @test round(vals[1], digits = 6) == 0.046933
-            @test round(vals[2], digits = 6) == 0.000841
-            @test round(vals[3], digits = 6) == 0.248696
+            @test round(vals[1], digits=6) == 0.046933
+            @test round(vals[2], digits=6) == 0.000841
+            @test round(vals[3], digits=6) == 0.248696
 
             # Probe (up and down are feasible)
             probe_up, probe_down = BB.probe(mip, mip.int_vars[1], 0.5, 0.0, 1.0, 1_000_000)
-            @test round(probe_down, digits = 6) == 62.690000
-            @test round(probe_up, digits = 6) == 62.714100
+            @test round(probe_down, digits=6) == 62.690000
+            @test round(probe_up, digits=6) == 62.714100
 
             # Fix one variable to zero
             BB.set_bounds!(mip, mip.int_vars[1:1], [0.0], [0.0])
             status, obj = BB.solve_relaxation!(mip)
             @test status == :Optimal
-            @test round(obj, digits = 6) == 62.690000
+            @test round(obj, digits=6) == 62.690000
 
             # Fix one variable to one and another variable variable to zero
             BB.set_bounds!(mip, mip.int_vars[1:2], [1.0, 0.0], [1.0, 0.0])
             status, obj = BB.solve_relaxation!(mip)
             @test status == :Optimal
-            @test round(obj, digits = 6) == 62.714777
+            @test round(obj, digits=6) == 62.714777
 
             # Fix all binary variables to one, making problem infeasible
             N = length(mip.int_vars)
@@ -68,7 +71,7 @@ function bb_run(optimizer_name, optimizer; large = true)
             BB.set_bounds!(mip, mip.int_vars, zeros(N), ones(N))
             status, obj = BB.solve_relaxation!(mip)
             @test status == :Optimal
-            @test round(obj, digits = 6) == 62.637280
+            @test round(obj, digits=6) == 62.637280
         end
 
         @testset "varbranch" begin
@@ -82,8 +85,8 @@ function bb_run(optimizer_name, optimizer; large = true)
                     BB.StrongBranching(),
                     BB.ReliabilityBranching(),
                     BB.HybridBranching(),
-                    BB.StrongBranching(aggregation = :min),
-                    BB.ReliabilityBranching(aggregation = :min, collect = true),
+                    BB.StrongBranching(aggregation=:min),
+                    BB.ReliabilityBranching(aggregation=:min, collect=true),
                 ]
                     h5 = H5File("$FIXTURES/$instance.h5")
                     mip_lower_bound = h5.get_scalar("mip_lower_bound")
@@ -98,23 +101,23 @@ function bb_run(optimizer_name, optimizer; large = true)
                     @info optimizer_name, branch_rule, instance
                     @time BB.solve!(
                         mip,
-                        initial_primal_bound = mip_primal_bound,
-                        print_interval = 1,
-                        node_limit = 25,
-                        branch_rule = branch_rule,
+                        initial_primal_bound=mip_primal_bound,
+                        print_interval=1,
+                        node_limit=25,
+                        branch_rule=branch_rule,
                     )
                 end
             end
         end
 
         @testset "collect" begin
-            rule = BB.ReliabilityBranching(collect = true)
+            rule = BB.ReliabilityBranching(collect=true)
             BB.collect!(
                 optimizer,
                 "$FIXTURES/bell5.mps.gz",
-                node_limit = 100,
-                print_interval = 10,
-                branch_rule = rule,
+                node_limit=100,
+                print_interval=10,
+                branch_rule=rule,
             )
             n_sb = rule.stats.num_strong_branch_calls
             h5 = H5File("$FIXTURES/bell5.h5")
@@ -132,3 +135,67 @@ function test_bb()
     @time bb_run("HiGHS", optimizer_with_attributes(HiGHS.Optimizer))
     # @time bb_run("CPLEX", optimizer_with_attributes(CPLEX.Optimizer, "CPXPARAM_Threads" => 1))
 end
+
+function test_bb_replay()
+    rule_sb = BB.StrongBranching()
+    rule_rb = BB.ReliabilityBranching()
+    optimizer = optimizer_with_attributes(HiGHS.Optimizer)
+    filenames = [replace(f, ".h5" => "") for f in glob("test/fixtures/stab/*.h5")]
+    results_filename = "tmp.csv"
+
+    lk = ReentrantLock()
+    results = []
+
+    function push_result(r)
+        lock(lk) do
+            push!(results, r)
+            df = DataFrame()
+            for row in results
+                push!(df, row, cols=:union)
+            end
+            CSV.write(results_filename, df)
+        end
+    end
+
+    function solve(filename; replay=nothing, skip=false, rule)
+        has_replay = (replay !== nothing)
+        h5 = H5File("$filename.h5", "r")
+        mip_obj_bound = h5.get_scalar("mip_obj_bound")
+        @show filename
+        @show has_replay
+        h5.file.close()
+        mip = BB.init(optimizer)
+        BB.read!(mip, "$filename.mps.gz")
+        time_solve = @elapsed begin
+            pool, replay = BB.solve!(
+                mip,
+                initial_primal_bound=mip_obj_bound,
+                print_interval=100,
+                node_limit=1_000,
+                branch_rule=rule,
+                replay=replay,
+            )
+        end
+        if !skip
+            push_result(
+                Dict(
+                    "Filename" => filename,
+                    "Replay?" => has_replay,
+                    "Solve time (s)" => time_solve,
+                    "Relative MIP gap (%)" => round(pool.gap * 100, digits=3)
+                )
+            )
+        end
+        return replay
+    end
+
+    # Solve reference instance
+    replay = solve(filenames[1], skip=true, rule=rule_sb)
+
+    # Solve perturbations
+    for i in 2:6
+        solve(filenames[i], rule=rule_rb, replay=nothing)
+        solve(filenames[i], rule=rule_rb, replay=deepcopy(replay))
+    end
+    return
+end