Add tests from Quantum Bridge Analytics

Project.toml

@@ -13,8 +13,9 @@ julia = "1.6"
 
 [extras]
 Anneal = "e4d9eb7f-b088-426e-aeb5-1c0dae3d8abb"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Anneal", "Test", "TOML"]
+test = ["Anneal", "JuMP", "Test", "TOML"]

src/MOI_wrapper.jl

@@ -25,9 +25,9 @@ function MOI.empty!(model::VirtualQUBOModel)
     isnothing(model.optimizer) || MOI.empty!(model.optimizer)
 
     # -*- PBFs -*-
-    empty!(model.ℍ)
-    empty!(model.ℍ₀)
-    empty!(model.ℍᵢ)
+    empty!(model.H)
+    empty!(model.H₀)
+    empty!(model.Hᵢ)
 
     # -*- MathOptInterface -*-
     empty!(model.moi)

src/ToQUBO.jl

@@ -47,4 +47,7 @@ include("qubo.jl")
 # -*- MOI Wrapper -*-
 include("MOI_wrapper.jl")
 
+
+include("analysis.jl")
+
 end # module

src/analysis.jl

@@ -0,0 +1,35 @@
+module Analysis
+    using ToQUBO: VirtualQUBOModel, PBO
+
+    function objective_function(model::VirtualQUBOModel)
+        model.H
+    end
+
+    function qubo_normal_form(model::VirtualQUBOModel)
+        PBO.qubo_normal_form(Array, objective_function(model))
+    end
+
+    function virtual_qubo_model(model)
+        if model isa VirtualQUBOModel
+            return model
+        end
+
+        if hasfield(typeof(model), :model)
+            inner = virtual_qubo_model(model.model)
+
+            if !isnothing(inner)
+                return inner
+            end
+        end
+
+        if hasfield(typeof(model), :optimizer)
+            inner = virtual_qubo_model(model.optimizer)
+
+            if !isnothing(inner)
+                return inner
+            end
+        end
+
+        return nothing
+    end
+end

src/model.jl

@@ -108,9 +108,9 @@ mutable struct VirtualQUBOModel{T} <: AbstractVirtualModel{T}
     # - Underlying Optimizer
     optimizer::Union{Nothing, MOI.AbstractOptimizer}
 
-    ℍ::ℱ{T}          # Total Energy
-    ℍ₀::ℱ{T}         # Objective
-    ℍᵢ::Vector{ℱ{T}} # Constraints
+    H::ℱ{T}          # Total Energy
+    H₀::ℱ{T}         # Objective
+    Hᵢ::Vector{ℱ{T}} # Constraints
 
     # -*- MathOptInterface -*-
     moi::VirtualQUBOModelMOI{T}

src/pbo.jl

@@ -612,9 +612,6 @@ function qubo_normal_form(::Type{<: AbstractDict}, f::PBF{S, T}) where {S, T}
 end
 
 function qubo_normal_form(::Type{<: AbstractArray}, f::PBF{S, T}) where {S, T}
-    # -* Constants *-
-    𝟐 = convert(T, 2)
-
     # -* QUBO *-
     x = varmap(f)
     n = length(x)
@@ -631,8 +628,8 @@ function qubo_normal_form(::Type{<: AbstractArray}, f::PBF{S, T}) where {S, T}
             Q[i, i] += a
         elseif k == 2
             i, j = η
-            Q[i, j] += a / 𝟐
-            Q[j, i] += a / 𝟐
+            Q[i, j] += a / 2
+            Q[j, i] += a / 2
         else
             error(DomainError, ": Can't convert Pseudo-boolean function with degree greater than 2 to QUBO format.\nTry using 'quadratize' before conversion.")
         end
@@ -646,14 +643,9 @@ function qubo_normal_form(f::PBF{S, T}) where {S, T}
     return qubo_normal_form(Dict, f)
 end
 
-function ising_normal_form(::Type{<: AbstractDict}, f::PBF{S, T}) where {S, T}
-    # -* Constants *-
-    𝟎 = zero(T)
-    𝟐 = convert(T, 2)
-    𝟒 = convert(T, 4)
-
+function ising_normal_form(::Type{<:AbstractDict}, f::PBF{S, T}) where {S, T}
     # -* QUBO *-
-    x, Q, c = qubo(Dict, f)
+    x, Q, c = qubo_normal_form(Dict, f)
 
     # -* Ising *-
     h = Dict{Int, T}()
@@ -663,18 +655,18 @@ function ising_normal_form(::Type{<: AbstractDict}, f::PBF{S, T}) where {S, T}
         i, j = ω
 
         if i == j
-            α = a / 𝟐
+            α = a / 2
 
-            h[i] = get(h, i, 𝟎) + α
+            h[i] = get(h, i, 0) + α
 
             c += α
         else
-            β = a / 𝟒
+            β = a / 4
 
             J[i, j] = β
 
-            h[i] = get(h, i, 𝟎) + β
-            h[j] = get(h, j, 𝟎) + β
+            h[i] = get(h, i, 0) + β
+            h[j] = get(h, j, 0) + β
 
             c += β
         end
@@ -684,12 +676,8 @@ function ising_normal_form(::Type{<: AbstractDict}, f::PBF{S, T}) where {S, T}
 end
 
 function ising_normal_form(::Type{<:AbstractArray}, f::PBF{S, T}) where {S, T}
-    # -* Constants *-
-    𝟐 = convert(T, 2)
-    𝟒 = convert(T, 4)
-
     # -* QUBO *-
-    x, Q, c = qubo(Dict, f)
+    x, Q, c = qubo_normal_form(Dict, f)
 
     # -* Ising *-
     n = length(x)
@@ -700,13 +688,13 @@ function ising_normal_form(::Type{<:AbstractArray}, f::PBF{S, T}) where {S, T}
         i, j = ω
 
         if i == j
-            α = a / 𝟐
+            α = a / 2
 
             h[i] += α
 
             c += α
         else
-            β = a / 𝟒
+            β = a / 4
 
             J[i, j] += β
 
@@ -775,7 +763,7 @@ end
 Defines a new quadratization technique.
 """
 macro quadratization(name, nsv, nst)
-    return :(
+    quote
         struct $(esc(name)) <: QuadratizationType end;
 
         function nsv(::Type{$(esc(name))}, k::Int)
@@ -785,7 +773,7 @@ macro quadratization(name, nsv, nst)
         function nst(::Type{$(esc(name))}, k::Int)
             return $(esc(nst))
         end;
-    )
+    end
 end
 
 @doc raw"""

src/qubo.jl

@@ -94,19 +94,19 @@ function toqubo!(model::VirtualQUBOModel{T}) where {T}
     # -*- :: Objective Function Assembly :: -*-
     ε = convert(T, 1.0) # TODO: This should be made a parameter too?
 
-    ρᵢ = δ(model.ℍ₀) ./ ϵ.(model.ℍᵢ; tol=model.settings.Tol) .+ ε
+    ρᵢ = δ(model.H₀) ./ ϵ.(model.Hᵢ; tol=model.settings.Tol) .+ ε
 
     if MOI.get(model, MOI.ObjectiveSense()) === MOI.MAX_SENSE
         ρᵢ *= -1.0
     end
 
-    model.ℍ = model.ℍ₀ + sum(ρᵢ .* model.ℍᵢ; init=zero(T))
+    model.H = model.H₀ + sum(ρᵢ .* model.Hᵢ; init=zero(T))
 
     Q = SQT{T}[]
     a = SAT{T}[]
     b = zero(T)
 
-    for (ω, c) in model.ℍ
+    for (ω, c) in model.H
         if length(ω) == 0
             b += c
         elseif length(ω) == 1
@@ -155,10 +155,10 @@ function toqubo_variables!(model::VirtualQUBOModel{T}) where {T}
     # ::: Variable Analysis :::
 
     # Set of all source variables
-    Ω = Set{VI}(MOI.get(model, MOI.ListOfVariableIndices()))
+    Ω = Vector{VI}(MOI.get(model, MOI.ListOfVariableIndices()))
 
     # Variable Sets and Bounds (Boolean, Integer, Real)
-    𝔹 = Set{VI}()
+    𝔹 = Vector{VI}()
     ℤ = Dict{VI, Tuple{Union{T, Nothing}, Union{T, Nothing}}}()
     ℝ = Dict{VI, Tuple{Union{T, Nothing}, Union{T, Nothing}}}()
 
@@ -280,7 +280,7 @@ function toqubo_objective!(model::VirtualQUBOModel{T}, F::Type{<:VI}) where {T}
     xᵢ = MOI.get(model, MOI.ObjectiveFunction{F}())
 
     for (ωᵢ, cᵢ) ∈ model.source[xᵢ]
-        model.ℍ₀[ωᵢ] += cᵢ
+        model.H₀[ωᵢ] += cᵢ
     end
 end
 
@@ -293,12 +293,12 @@ function toqubo_objective!(model::VirtualQUBOModel{T}, F::Type{<:SAF{T}}) where
         xᵢ = aᵢ.variable
 
         for (ωᵢ, dᵢ) ∈ model.source[xᵢ]
-            model.ℍ₀[ωᵢ] += cᵢ * dᵢ
+            model.H₀[ωᵢ] += cᵢ * dᵢ
         end
     end
 
     # -*- Constant -*-
-    model.ℍ₀ += f.constant
+    model.H₀ += f.constant
 end
 
 function toqubo_objective!(model::VirtualQUBOModel{T}, F::Type{<:SQF{T}}) where {T}
@@ -311,8 +311,15 @@ function toqubo_objective!(model::VirtualQUBOModel{T}, F::Type{<:SQF{T}}) where
         xᵢ = Qᵢ.variable_1
         xⱼ = Qᵢ.variable_2
 
+        # MOI convetion is to write ScalarQuadraticFunction as
+        #     ½ x' Q x + a x + b
+        # ∴ every coefficient in the main diagonal is doubled
+        if xᵢ === xⱼ
+            cᵢ /= 2
+        end
+
         for (ωᵢ, dᵢ) ∈ model.source[xᵢ], (ωⱼ, dⱼ) ∈ model.source[xⱼ]
-            model.ℍ₀[ωᵢ × ωⱼ] += cᵢ * dᵢ * dⱼ
+            model.H₀[ωᵢ × ωⱼ] += cᵢ * dᵢ * dⱼ
         end
     end
 
@@ -321,12 +328,12 @@ function toqubo_objective!(model::VirtualQUBOModel{T}, F::Type{<:SQF{T}}) where
         xᵢ = aᵢ.variable
 
         for (ωᵢ, dᵢ) in model.source[xᵢ]
-            model.ℍ₀[ωᵢ] += cᵢ * dᵢ
+            model.H₀[ωᵢ] += cᵢ * dᵢ
         end
     end
 
     # -*- Constant -*-
-    model.ℍ₀ += f.constant
+    model.H₀ += f.constant
 end
 
 @doc raw"""
@@ -356,7 +363,7 @@ function toqubo_constraint!(model::VirtualQUBOModel{T}, F::Type{<:SAF{T}}, S::Ty
         gᵢ = PBO.discretize((gᵢ - bᵢ) ^ 2; tol=model.settings.Tol)
         hᵢ = PBO.quadratize(gᵢ; slack = slack_factory(model))
 
-        push!(model.ℍᵢ, hᵢ)
+        push!(model.Hᵢ, hᵢ)
     end
 
     nothing
@@ -389,7 +396,7 @@ function toqubo_constraint!(model::VirtualQUBOModel{T}, F::Type{<:SAF{T}}, S::Ty
         sᵢ = ℱ{T}(collect(slackℤ!(model; α=αᵢ, β=βᵢ, name=:s)))
         hᵢ = PBO.quadratize((gᵢ + sᵢ) ^ 2;slack = slack_factory(model))
 
-        push!(model.ℍᵢ, hᵢ)
+        push!(model.Hᵢ, hᵢ)
     end
 
     nothing
@@ -408,6 +415,10 @@ function toqubo_constraint!(model::VirtualQUBOModel{T}, F::Type{<:SQF{T}}, S::Ty
             xⱼ = Qⱼ.variable_1
             yⱼ = Qⱼ.variable_2
 
+            if xᵢ === xⱼ
+                cᵢ /= 2
+            end
+
             for (ωⱼ, dⱼ) ∈ model.source[xⱼ], (ηⱼ, eⱼ) ∈ model.source[yⱼ]
                 gᵢ[ωⱼ × ηⱼ] += cⱼ * dⱼ * eⱼ
             end
@@ -425,7 +436,7 @@ function toqubo_constraint!(model::VirtualQUBOModel{T}, F::Type{<:SQF{T}}, S::Ty
         gᵢ = PBO.discretize((gᵢ - bᵢ) ^ 2; tol=model.settings.Tol)
         hᵢ = PBO.quadratize(gᵢ; slack = slack_factory(model))
 
-        push!(model.ℍᵢ, hᵢ)
+        push!(model.Hᵢ, hᵢ)
     end
 
     nothing
@@ -444,6 +455,10 @@ function toqubo_constraint!(model::VirtualQUBOModel{T}, F::Type{<:SQF{T}}, S::Ty
             xⱼ = Qⱼ.variable_1
             yⱼ = Qⱼ.variable_2
 
+            if xᵢ === xⱼ
+                cᵢ /= 2
+            end
+
             for (ωⱼ, dⱼ) ∈ model.source[xⱼ], (ηⱼ, eⱼ) ∈ model.source[yⱼ]
                 gᵢ[ωⱼ × ηⱼ] += cⱼ * dⱼ * eⱼ
             end
@@ -467,7 +482,7 @@ function toqubo_constraint!(model::VirtualQUBOModel{T}, F::Type{<:SQF{T}}, S::Ty
         sᵢ = ℱ{T}(collect(slackℤ!(model; α=αᵢ, β=βᵢ, name=:s)))
         hᵢ = PBO.quadratize((gᵢ + sᵢ) ^ 2; slack = slack_factory(model))
 
-        push!(model.ℍᵢ, hᵢ)
+        push!(model.Hᵢ, hᵢ)
     end
 end
 

test/examples/qba2.jl

@@ -0,0 +1,34 @@
+@testset "Illustrative Example" begin
+
+Q̄ = [-5  2  4  0
+      2 -3  1  0
+      4  1 -8  5 
+      0  0  5 -6 ]
+
+c̄ = 0
+x̄ = [1, 0, 0, 1]
+ȳ = -11
+
+model = Model(() -> ToQUBO.Optimizer(ExactSampler.Optimizer))
+
+@variable(model, x[i = 1:4], Bin)
+@objective(model, Min, -5x[1] - 3x[2] - 8x[3] - 6x[4] + 4x[1] * x[2] + 8x[1] * x[3] + 2x[2] * x[3] + 10x[3] * x[4])
+
+optimize!(model)
+
+vqm = unsafe_backend(model)
+
+_, Q, c = ToQUBO.Analysis.qubo_normal_form(vqm)
+
+x̂ = value.(x)
+ŷ = objective_value(model)
+
+# :: Reformulation ::
+@test c == c̄
+@test Q == Q̄
+
+# :: Solution ::
+@test x̂ == x̄
+@test ŷ == ȳ
+
+end