Benchmarking¶
David E. Bernal NeiraDavidson School of Chemical Engineering, Purdue University
Universities Space Research Association
NASA QuAIL
Pedro Maciel Xavier
Davidson School of Chemical Engineering, Purdue University
Environment Setup¶
If you open this notebook in Google Colab, make sure the runtime is set to Julia before running the next cell.
The setup cell will clone SECQUOIA/QuIP into the Colab workspace when needed, activate notebooks_jl/Project.toml, and install the Julia packages used in this notebook.
import Pkg
IN_COLAB = haskey(ENV, "COLAB_RELEASE_TAG") || haskey(ENV, "COLAB_JUPYTER_IP")
function detect_quip_repo_dir()
candidates = (pwd(), normpath(pwd(), ".."))
for candidate in candidates
if isfile(joinpath(candidate, "notebooks_jl", "Project.toml"))
return candidate
end
end
error("Could not locate the QuIP repository root from $(pwd()).")
end
QUIP_REPO_DIR = if IN_COLAB
repo_dir = get(ENV, "QUIP_REPO_DIR", joinpath(pwd(), "QuIP"))
if !isdir(repo_dir)
run(`git clone --depth 1 https://github.com/SECQUOIA/QuIP.git $repo_dir`)
end
repo_dir
else
detect_quip_repo_dir()
end
JULIA_NOTEBOOKS_DIR = joinpath(QUIP_REPO_DIR, "notebooks_jl")
Pkg.activate(JULIA_NOTEBOOKS_DIR)
Pkg.instantiate(; io = devnull)
cd(JULIA_NOTEBOOKS_DIR)
About this Notebook¶
With the new availability of unconventional hardware, novel algorithms, and increasingly optimized software to address optimization problems; the first question that arises is, which one is better? We will call the combination of hardware, algorithm, and software within a solution method as a solver. This question is also relevant when evaluating a single solver, given that usually they rely on hyperparameters, for which the quesiton now becomes into, which is the best parameter setting for a given solver? These questions obviously depend on the problem that one is trying to solve. The solution of the problem also depends on the budget of resources that one has available.
In the case that the available resources are relatively “unlimited” and that the problem to solve is known, one could exhaustively try all the parameter settings within a delimited range for that instance and choose which one is the best. This case is idealistic, in the sense that usually one does not know a-priori which problem is there to solve (and if while testing all the parameters you solve it, what would be the point of identifying the best parameters?), and that there exists limitations in terms of resources, e.g., time, memory, or energy, when trying to address these problems. A case closer to reality is where you have the chance of solving problems that look similar to the one that you are interested in solving later, either because you have previously generated problems or you have identified a feature that characterizes your problem of interest and can generate random instances, which we will call as a family of instances. Then, you can use a larger amount of resources to solve that family of problems “off-line”, meaning that you spend extra resources to address the problems in your family of instances although it is unrelated to the actual application. Finally, you would like to use the results that you found off-line as a guidance to solve your unknown problem more efficiently.
Example¶
For illustration purposes, we will use an example that you are already familiar with, which is an Ising model. As a solver, we will use the PySA simulated annealing code.
Ising model¶
This notebook will explain the basics of the Ising model. In order to implement the different Ising Models we will use JuMP and QUBO.jl, for defining the Ising model and solving it with simulated annealing.
Problem statement¶
We pose the Ising problem as the following optimization problem:
where we optimize over spins , on a constrained graph , where the quadratic coefficients are and the linear coefficients are . We also include an arbitrary offset of the Ising model .
using JuMP
using QUBO
using LinearAlgebrah = [
145.0,
122.0,
122.0,
266.0,
266.0,
266.0,
242.5,
266.0,
386.5,
387.0,
386.5,
]
J = [
0 0 0 24 24 24 0 24 24 24 24
0 0 0 24 0 24 24 0 24 24 24
0 0 0 0 24 0 24 24 24 24 24
0 0 0 0 24 48 24 24 48 48 48
0 0 0 0 0 24 24 48 48 48 48
0 0 0 0 0 0 24 24 48 48 48
0 0 0 0 0 0 0 24 48 48 48
0 0 0 0 0 0 0 0 48 48 48
0 0 0 0 0 0 0 0 0 72 72
0 0 0 0 0 0 0 0 0 0 72
0 0 0 0 0 0 0 0 0 0 0
]
β = 1319.5
ising_model = Model()
@variable(ising_model, s[1:11], Spin)
@objective(ising_model, Min, s' * J * s + h' * s + β)
println(ising_model)Min 24 s[4]*s[1] + 24 s[4]*s[2] + 24 s[5]*s[1] + 24 s[5]*s[3] + 24 s[5]*s[4] + 24 s[6]*s[1] + 24 s[6]*s[2] + 48 s[6]*s[4] + 24 s[6]*s[5] + 24 s[7]*s[2] + 24 s[7]*s[3] + 24 s[7]*s[4] + 24 s[7]*s[5] + 24 s[7]*s[6] + 24 s[8]*s[1] + 24 s[8]*s[3] + 24 s[8]*s[4] + 48 s[8]*s[5] + 24 s[8]*s[6] + 24 s[8]*s[7] + 24 s[9]*s[1] + 24 s[9]*s[2] + 24 s[9]*s[3] + 48 s[9]*s[4] + 48 s[9]*s[5] + 48 s[9]*s[6] + 48 s[9]*s[7] + 48 s[9]*s[8] + 24 s[10]*s[1] + 24 s[10]*s[2] + 24 s[10]*s[3] + 48 s[10]*s[4] + 48 s[10]*s[5] + 48 s[10]*s[6] + 48 s[10]*s[7] + 48 s[10]*s[8] + 72 s[10]*s[9] + 24 s[11]*s[1] + 24 s[11]*s[2] + 24 s[11]*s[3] + 48 s[11]*s[4] + 48 s[11]*s[5] + 48 s[11]*s[6] + 48 s[11]*s[7] + 48 s[11]*s[8] + 72 s[11]*s[9] + 72 s[11]*s[10] + 145 s[1] + 122 s[2] + 122 s[3] + 266 s[4] + 266 s[5] + 266 s[6] + 242.5 s[7] + 266 s[8] + 386.5 s[9] + 387 s[10] + 386.5 s[11] + 1319.5
Subject to
s[1] spin
s[2] spin
s[3] spin
s[4] spin
s[5] spin
s[6] spin
s[7] spin
s[8] spin
s[9] spin
s[10] spin
s[11] spin
We can visualize the graph that defines this instance using the matrix as the adjacency matrix of a graph.
using Plots
# Make plots look professional
Plots.default(;
fontfamily = "Computer Modern",
plot_titlefontsize = 16,
titlefontsize = 14,
guidefontsize = 12,
legendfontsize = 10,
tickfontsize = 10,
)# plot(QUBOTools.SystemLayoutPlot(J + diagm(h)))
plot(QUBOTools.backend(unsafe_backend(ising_model)))Cannot convert QUBOTools.Model{MathOptInterface.VariableIndex, Float64, Int64} to series data for plotting
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] _prepare_series_data(x::QUBOTools.Model{MathOptInterface.VariableIndex, Float64, Int64})
@ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/series.jl:8
[3] _series_data_vector(x::QUBOTools.Model{MathOptInterface.VariableIndex, Float64, Int64}, plotattributes::Dict{Symbol, Any})
@ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/series.jl:36
[4] macro expansion
@ ~/.julia/packages/RecipesPipeline/BGM3l/src/series.jl:129 [inlined]
[5] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, #unused#::Type{RecipesPipeline.SliceIt}, x::Any, y::Any, z::Any)
@ RecipesPipeline ~/.julia/packages/RecipesBase/BRe07/src/RecipesBase.jl:300
[6] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
@ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/user_recipe.jl:38
[7] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
@ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/RecipesPipeline.jl:72
[8] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
@ Plots ~/.julia/packages/Plots/sxUvK/src/plot.jl:223
[9] #plot#188
@ ~/.julia/packages/Plots/sxUvK/src/plot.jl:102 [inlined]
[10] plot(args::Any)
@ Plots ~/.julia/packages/Plots/sxUvK/src/plot.jl:93
[11] top-level scope
@ /mnt/c/Users/pedroxavier/QUBONotebooks/notebooks/5-Benchmarking.ipynb:2Since the problem is relatively small (11 variables, combinations), we can afford to enumerate all the solutions.
set_optimizer(ising_model, QUBO.ExactSampler.Optimizer)
optimize!(ising_model)
enumeration_time = solve_time(ising_model)
@info "Enumeration took $(enumeration_time) seconds"┌ Info: Enumeration took 0.000969682 seconds
└ @ Main /mnt/c/Users/pedroxavier/QUBONotebooks/notebooks/5-Benchmarking.ipynb:7
using DWaveNealset_optimizer(ising_model, DWaveNeal.Optimizer)
set_optimizer_attribute(ising_model, "num_reads", 1_000)
optimize!(ising_model)
ising_s = round.(Int, value.(s))
# Display solution of the problem
println(solution_summary(ising_model))
println("* s = $ising_s")* Solver : D-Wave Neal Simulated Annealing Sampler
* Status
Result count : 12
Termination status : LOCALLY_SOLVED
Message from the solver:
""
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : NO_SOLUTION
Objective value : 5.00000e+00
* Work counters
Solve time (sec) : 2.70403e-01
* s = [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 1]
# plot(QUBOTools.EnergyFrequencyPlot(ising_model))
plot(QUBOTools.sampleset(unsafe_backend(ising_model)))We are going to use the default limits of temperature given by the simulating annealing code. These are defined using the minimum and maximum nonzero coefficients in the Ising model. Then the range for beta is defined as
where
Hot temperature: We want to scale hot_beta so that for the most unlikely qubit flip, we get at least 50% chance of flipping. (This means all other qubits will have > 50% chance of flipping initially). Most unlikely flip is when we go from a very low energy state to a high energy state, thus we calculate hot_beta based on max_delta_energy.
Cold temperature: Towards the end of the annealing schedule, we want to minimize the chance of flipping. Don’t want to be stuck between small energy tweaks. Hence, set cold_beta so that at minimum energy change, the chance of flipping is set to 1%.
By default, the schedule also follows a geometric series.
function geomspace(a, b; length = 100)
return exp10.(range(log10(a), log10(b); length))
endgeomspace (generic function with 1 method)# data = QUBOTools.metadata(QUBOTools.solution(qubo_model))
# info = data["info"]
β₀, β₁ = (2.0, 100.0) # info["beta_range"](2.0, 100.0)function plot_schedule(β₀, β₁; length = 1_000)
β = geomspace(β₀, β₁; length=1_000)
plt = plot(;
plot_title = "Default Geometric temperature schedule",
xlabel = "Sweeps",
ylabel = raw"$ \beta $ = Inverse temperature",
)
plot!(
plt, β;
color = cgrad(:redblue; rev=true),
linestyle = :dot,
linewidth = 3,
line_z = β,
legend = nothing,
yticks = ([β₀, β₁], [raw"$ \beta_{0} $", raw"$ \beta_{1} $"]),
)
return plt
end
plot_schedule(β₀, β₁)sweeps = [1:250; 260:10:1_000]
schedules = ["geometric", "linear"]
s = 0.99 # Success probability
λ⃰ = 5.0 # Optimal Energy \lambda[tab]\asteraccent[tab]
results = Dict{Symbol,Any}(
:p => Dict{String,Any}(), # Success probability
:t => Dict{String,Any}(), # Time elapsed
:ttt => Dict{String,Any}(), # Time-to-target
)
set_attribute(ising_model, "num_reads", 1_000)
for schedule in schedules
results[:p][schedule] = Float64[]
results[:t][schedule] = Float64[]
results[:ttt][schedule] = Float64[]
set_attribute(ising_model, "beta_schedule_type", schedule)
for sweep in sweeps
set_attribute(ising_model, "num_sweeps", sweep)
optimize!(ising_model)
sol = QUBOTools.sampleset(unsafe_backend(ising_model))
p = QUBOTools.success_rate(sol, λ⃰)
t = QUBOTools.effective_time(sol) # TODO: Should be 'total_time'?
ttt = QUBOTools.tts(sol, λ⃰, s)
push!(results[:p][schedule], p)
push!(results[:t][schedule], t)
push!(results[:ttt][schedule], ttt)
end
end
resultsDict{Symbol, Any} with 3 entries:
:p => Dict{String, Any}("geometric"=>[0.64, 0.507, 0.426, 0.397, 0.335, 0.3…
:ttt => Dict{String, Any}("geometric"=>[0.0126028, 0.019915, 0.0227079, 0.026…
:t => Dict{String, Any}("geometric"=>[0.00279592, 0.00305847, 0.00273731, 0…function plot_schedule_grid(sweeps, schedules, results)
plt = plot(;
layout = @layout([a{0.475h}; b{0.475h}; leg{0.05h}]),
size = (600, 600),
plot_title = raw"""
Simulated annealing expected runtime of easy
$ n = 11 $ Ising model with varying schedule and sweeps
""",
plot_titlevspan = 0.1,
)
# Elapsed time
plot!(
plt;
subplot = 1,
ylabel = "Time [s]",
yscale = :log10,
legend = false,
)
for schedule in schedules
plot!(
plt, sweeps, results[:t][schedule];
subplot = 1,
linestyle = :solid,
label = schedule,
linewidth = 3,
)
end
hline!(
plt, [results[:t]["geometric"][end]];
subplot = 1,
linestyle = :dash,
label = "default",
color = :gray,
linewidth = 1,
)
hline!(
plt, [enumeration_time];
subplot = 1,
linestyle = :dash,
label = "enumeration",
color = :black,
linewidth = 1,
)
# Success probability
plot!(
plt;
subplot = 2,
xlabel = "Sweeps",
ylabel = "Success Probability [%]",
legend = false,
)
for schedule in schedules
plot!(
plt, sweeps, results[:p][schedule];
subplot = 2,
linestyle = :solid,
label = schedule,
ymirror = true,
linewidth = 3,
)
end
hline!(
plt, [results[:p]["geometric"][end]];
subplot = 2,
linestyle = :dash,
label = "default",
color = :gray,
linewidth = 1,
)
hline!(
plt, [1.0];
subplot = 2,
linestyle = :dash,
label = "enumeration",
color = :black,
linewidth = 1,
)
# Legend
plot!(
plt, (1:4)';
subplot = 3,
framestyle = nothing,
showaxis = false,
grid = false,
linestyle = [:solid :solid :dash :dash],
color = [1 2 :gray :black],
label = ["geometric" "linear" "default" "enumeration"],
linewidth = [3 3 1 1],
legend = :top,
legend_column = 4,
)
return plt
end
plot_schedule_grid(sweeps, schedules, results)These plots represent often contradictory metrics: on one hand you would like to obtain a large probability of finding a right solution (the definition of right comes from what you define as success). On the other hand, the time it takes to solve these cases should be as small as possible. This is why we are interested in a metric that combines both, and that is why we settle on the Time To Solution (TTS) which is defined as
where s is a success factor, usually takes as , and is the success probability, usually accounted as the observed success probability.
One usually reads this as the time to solution within probability.
function plot_ttt_grid(sweeps, schedules, results, k = 20)
plt = plot(;
layout = @layout([a{0.95h}; leg{0.05h}]),
size = (700, 600),
plot_title = raw"""
Simulated annealing expected runtime of easy
$ n = 11 $ Ising model with varying schedule and sweeps
""",
plot_titlevspan = 0.1,
plot_titlefont = Plots.font(16, "Computer Modern"),
titlefont = Plots.font(14, "Computer Modern"),
guidefont = Plots.font(12, "Computer Modern"),
legendfont = Plots.font(10, "Computer Modern"),
tickfont = Plots.font(10, "Computer Modern"),
inset_subplots = [(1, bbox(0.25,0.166,0.4,0.4, :center))],
)
# Elapsed time
plot!(
plt;
subplot = 1,
xlabel = "Sweeps",
ylabel = "Time to target with \$ $(100s) \\% \$ chance [s]",
legend = false,
legend_column = 3,
)
for schedule in schedules
plot!(
plt, sweeps[k+1:end], results[:ttt][schedule][k+1:end];
subplot = 1,
linestyle = :solid,
label = schedule,
yscale = :log10,
linewidth = 3,
)
end
# Value for the default solution
hline!(
plt, [results[:ttt]["geometric"][end]];
subplot = 1,
linestyle = :dash,
label = "default",
color = :gray,
linewidth = 1,
)
# Legend
plot!(
plt, (1:3)';
subplot = 2,
framestyle = nothing,
showaxis = false,
grid = false,
linestyle = [:solid :solid :dash],
linewidth = [3 3 1],
color = [1 2 :gray],
label = ["geometric" "linear" "default"],
legend = :top,
legend_column = 3,
)
plot!(
plt[3];
xlabel = "Sweeps",
ylabel = "TTT [s]",
yscale = :log10,
legend = false,
framestyle = :box,
)
for schedule in schedules
plot!(
plt[3], sweeps[1:k], results[:ttt][schedule][1:k];
linewidth = 2,
marker = :diamond,
markersize = 4,
)
end
hline!(
plt[3], [results[:ttt]["geometric"][end]];
linestyle = :dash,
color = :gray,
linewidth = 1,
)
return plt
end
plot_ttt_grid(sweeps, schedules, results)As you can notice, the default parameters given by D-Wave (number of sweeps = 1000 and a geometric update of ) are not optimal for our tiny example in terms of expected runtime. This is certainly a function of the problem, for such a small instance having two sweeps are more than enough and more sweeps are an overkill. This parameters choice might not generalize to any other problem, as seen below.
Example 2¶
Let’s define a larger model, with 100 variables and random weights, to see how this performance changes.
Assume that we are interested at the instance created with random weights .
using Random# Number of variables
n = 100
# Fixing the random seed to get the same result
Random.seed!(1)
# We only consider upper triangular matrix ignoring the diagonal
J = triu!(2rand(Float64, n, n) .- 1, 1)
h = 2rand(Float64, n) .- 1
;plot(QUBOTools.SystemLayoutPlot(J + diagm(h)))For a problem of this size we cannot do a complete enumeration () but we can randomly sample the distribution of energies to have a baseline for our later comparisons.
random_ising_model = Model()
@variable(random_ising_model, s[1:n], Spin)
@objective(random_ising_model, Min, s' * J * s + h' * s)
random_ising_modelA JuMP Model
Minimization problem with:
Variables: 100
Objective function type: QuadExpr
`VariableRef`-in-`Spin`: 100 constraints
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
Names registered in the model: susing Statisticsset_optimizer(random_ising_model, QUBO.RandomSampler.Optimizer)
set_attribute(random_ising_model, "num_reads", 1_000)
optimize!(random_ising_model)
random_sampling_time = solve_time(random_ising_model)
energies = [objective_value(random_ising_model; result = i) for i = 1:result_count(random_ising_model)]
random_mean_energy = mean(energies)
@info "Average random energy = $(random_mean_energy)"┌ Info: Average random energy = 1.100341856327635
└ @ Main /mnt/c/Users/pedroxavier/QUBONotebooks/notebooks/5-Benchmarking.ipynb:12
set_optimizer(random_ising_model, DWaveNeal.Optimizer)
set_attribute(random_ising_model, "num_reads", 1_000)
optimize!(random_ising_model)
@info "Simulated Annealing best energy = $(objective_value(random_ising_model))"┌ Info: Simulated Annealing best energy = -408.47197451530747
└ @ Main /mnt/c/Users/pedroxavier/QUBONotebooks/notebooks/5-Benchmarking.ipynb:6
Notice that the minimum energy coming from the random sampling and the one from the simulated annealing are very different. Moreover, the distributions that both lead to are extremely different too.
# data = QUBOTools.metadata(QUBOTools.solution(qubo_model))
# info = data["info"]
β₀, β₁ = (2.0, 100.0) # info["beta_range"]
plot_schedule(β₀, β₁)We can solve this problem using IP such that we have guarantees that it is solved to optimality (this might be a great quiz for future lectures), but in this case let us define the “success” as getting an objective certain percentage of the best found solution in all cases (which we see it might not be even found with the default parameters). To get a scaled version of this success equivalent for all instances, we will define this success with respect to the metric:
Where corresponds to the best found solution within our sampling, is the mean of the random sampling shown above, and corresponds to the best found solution to our problem during the exploration. Consider that this minimum might not be the global minimum. This metric is very informative given that the best performance you can have is 1, being at the minimum, and negative values would correspond to a method that at best behaves worse that the random sampling. Success now is counted as being within certain treshold of this value of 1. This new way of measuring each method is very similar to the approximation ratio of approximation algorithms, therefore we will use that terminology from now on.
Before figuring out if we have the right optimal parameters, we want to save some effort by loading previously computed results.
If you do not want to load the results that we are providing, feel free to change the overwrite_pickles variable, at the expense that it will take some time (around 3 minutes per instance) to run.
If you do not want to wait, drop the results.zip file in the folder that is about to be created.
pickle_path = joinpath(JULIA_NOTEBOOKS_DIR, "results")
if !isdir(pickle_path)
@warn "Results directory '$pickle_path' does not exist. We will create it."
mkpath(pickle_path)
end
If you have a results.zip archive with precomputed Julia samples, place it in that folder and the next cell will unpack it.
Otherwise the notebook will generate the benchmarking data locally when needed.
zip_name = joinpath(pickle_path, "results.zip")
overwrite_pickles = false
if isfile(zip_name)
run(`unzip -o $zip_name -d $pickle_path`)
@info("Results archive extracted to '$pickle_path'")
end
use_raw_data = any(endswith(name, ".p") for name in readdir(pickle_path))
Now either we have the pickled file or not, let us compute the statistics we are looking for.
s = 0.99 # This is the success probability for the ttt calculation
threshold = 5.0 # This is a percentual threshold of what the minimum energy should be
sweeps = [1:250;260:10:1_000]
schedules = ["geometric", "linear"]
total_reads = 1_000
default_sweeps = 1_000
n_boot = 500
ci = 68 # Confidence interval for bootstrapping
default_boots = default_sweeps
boots = [1, 10, default_boots]
min_energy = -239.7094652034834
instance = 42
results_name = "results_$(instance).json"
results_name = joinpath(pickle_path, results_name)
results = Dict{Symbol,Any}(
:p => Dict{Int,Any}(),
:min_energy => Dict{Int,Any}(),
:random_energy => Dict{Int,Any}(),
:ttt => Dict{Int,Any}(),
:tttci => Dict{Int,Any}(),
:t => Dict{Int,Any}(),
:best => Dict{Int,Any}(),
:bestci => Dict{Int,Any}(),
)
# If you wanto to use the raw data and process it here
if use_raw_data || !isfile(results_name)
# If you want to generate the data or load it here
overwrite_pickles = false
for boot in boots
results[:p][boot] = Dict{X,Any}()
results[:ttt][boot] = Dict{X,Any}()
results[:tttci][boot] = Dict{X,Any}()
results[:best][boot] = Dict{X,Any}()
results[:bestci][boot] = Dict{X,Any}()
end
for schedule in schedules
probs = Dict{Int,Any}(k => [] for k in boots)
time_to_sol = Dict{Int,Any}(k => [] for k in boots)
prob_np = Dict{Int,Any}(k => [] for k in boots)
tttcs = Dict{Int,Any}(k => [] for k in boots)
times = []
b = Dict{Int,Any}(k => [] for k in boots)
bnp = Dict{Int,Any}(k => [] for k in boots)
bcs = Dict{Int,Any}(k => [] for k in boots)
for sweep in sweeps
# Gather instance names
pickle_name = joinpath(pickle_path, "$instance)_$(schedule)_$(sweep).p")
# If the instance data exists, load the data
if isfile(pickle_name) && !overwrite_pickles
# print(pickle_name)
sol = QUBOTools.read_solution(pickle_name)
time_s = QUBOTools.effective_time(sol)
else # If it does not exist, generate the data
set_attribute(random_ising_model, "num_reads", total_reads)
set_attribute(random_ising_model, "num_sweeps", sweep)
set_attribute(random_ising_model, "beta_schedule_type", schedule)
optimize!(random_ising_model)
sol = QUBOTools.solution(random_ising_model)
QUBOTools.write_solution(pickle_name, sol)
end
# Compute statistics
energies = QUBOTools.value.(sol)
occurrences = QUBOTools.reads.(sol)
total_counts = sum(occurrences)
push!(times, time_s)
if min(energies) < min_energy
min_energy = min(energies)
@info("A better solution of '$(min_energy)' was found for sweep '$sweep'")
end
# success = min_energy*(1.0 + threshold/100.0)**sign(min_energy)
success = random_energy - (random_energy - min_energy)*(1.0 - threshold/100.0)
# Best of boot samples es computed via n_boot bootstrappings
boot_dist = Dict{Int,Any}()
pr_dist = Dict{Int,Any}()
cilo = Dict{Int,Any}()
ciup = Dict{Int,Any}()
pr = Dict{Int,Any}()
pr_cilo = Dict{Int,Any}()
pr_ciup = Dict{Int,Any}()
for boot in boots
boot_dist[boot] = []
pr_dist[boot] = []
for i in 1:n_boot
resampler = rand(0:total_reads, boot)
sample_boot = energies.take(resampler, axis=0)
# Compute the best along that axis
push!(boot_dist[boot], min(sample_boot))
occurences = occurrences[resampler]
counts = Dict{Float64,Int}()
for (index, energy) in enumerate(sample_boot)
if energy in keys(counts)
counts[energy] += occurences[index]
else
counts[energy] = occurences[index]
end
end
push!(
pr_dist[boot],
sum(counts[key] for key in keys(counts) if key < success) / boot
)
end
push!(b[boot], mean(boot_dist[boot]))
# Confidence intervals from bootstrapping the best out of boot
bnp[boot] = boot_dist[boot]
cilo[boot] = map(x -> percentile(x, 50 - ci / 2), bnp[boot])
ciup[boot] = map(x -> percentile(x, 50 + ci / 2), bnp[boot])
push!(bcs[boot], (cilo[boot], ciup[boot]))
# Confidence intervals from bootstrapping the ttt of boot
prob_np[boot] = pr_dist[boot]
pr[boot] = mean(prob_np[boot])
push!(probs[boot], pr[boot])
if !all(x -> x > 0, prob_np[boot])
push!(time_to_sol[boot], Inf)
push!(tttcs[boot], (Inf, Inf))
else
pr_cilo[boot] = map(x -> percentile(x, 50 - ci / 2), prob_np[boot])
pr_ciup[boot] = map(x -> percentile(x, 50 + ci / 2), prob_np[boot])
push!(time_to_sol[boot], time_s * log10(1-s) / log10(1 - pr[boot] + 1E-9))
push!(
tttcs[boot],
(
time_s*log10(1-s)/log10(1 - pr_cilo[boot]),
time_s*log10(1-s)/log10(1 - pr_ciup[boot]+1E-9)
)
)
end
end
end
results[:t][schedule] = times
results[:min_energy][schedule] = min_energy
results[:random_energy][schedule] = random_energy
for boot in boots
results[:p][boot][schedule] = probs[boot]
results[:ttt][boot][schedule] = time_to_sol[boot]
results[:tttci][boot][schedule] = tttcs[boot]
results[:best][boot][schedule] = [
(random_energy - energy) / (random_energy - min_energy)
for energy in b[boot]
]
results[:bestci][boot][schedule] = [
tuple((random_energy - element) / (random_energy - min_energy) for element in energy)
for energy in bcs[boot]
]
end
end
# Save results file in case that we are interested in reusing them
JSON.print(results_name, results)
else # Just reload processed datafile
results = JSON.parsefile(results_name)
end
After gathering all the results, we would like to see the progress of the approximation ratio with respect to the increasing number of sweeps. To account for the stochasticity of this method, we are bootstrapping all of our results with different values of the bootstrapping sample, and each confidence interval corresponds to a standard deviation away from the mean.
function plot_progress(sweeps, boots, schedules, results)
plt = plot(;
plot_title = """
Simulated annealing approximation ratio of Ising 42 N=100
with varying schedule, $(n_boot) bootstrap re-samples, and sweeps
""",
xscale = :log10,
# xlims = (1, 200),
ylims = (0.8, 1.01),
xlabel = "Sweeps",
ylabel = raw"Approximation ratio = \
$\frac{\textrm{best found} - \textrm{random sample}}\
{\textrm{min energy} - \textrm{random sample}}$",
legend = :outertop,
)
for boot in boots
for schedule in schedules
bestnp = transpose(stack(results[:bestci][boot][schedule], dims=1))
plot!(
plt,
sweeps, bestnp[0];
fillrange = bestnp[1],
fillalpha = 0.25,
)
plot!(
plt,
sweeps, results[:best][boot][schedule];
label = "$schedule, $boot reads",
)
end
end
end
plot_progress (generic function with 1 method)Ensemble Analysis Note¶
The original lecture material included a later benchmarking comparison block that still existed only in Python. That mixed-language tail has been removed from this Julia notebook so the Colab and Jupyter workflow stays executable from top to bottom. A Julia-native ensemble analysis section can be added back in a later cleanup pass.