diff --git a/Project.toml b/Project.toml index 19b55f53..530e9dfc 100644 --- a/Project.toml +++ b/Project.toml @@ -1,17 +1,21 @@ -name = "ToQUBO" -uuid = "9a412ddf-83fa-43b6-9748-7843c851aa65" -authors = ["pedromxavier ", "joaquimg ", "AndradeTiago ", "David Bernal "] +name = "ToQUBO" +uuid = "9a412ddf-83fa-43b6-9748-7843c851aa65" +authors = [ + "pedromxavier ", + "pedroripper ", + "joaquimg ", + "AndradeTiago ", + "bernalde ", +] version = "0.1.4" [deps] -MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" -MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" -QUBOTools = "60eb5b62-0a39-4ddc-84c5-97d2adff9319" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +QUBOTools = "60eb5b62-0a39-4ddc-84c5-97d2adff9319" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" [compat] -MathOptInterface = "1" -MutableArithmetics = "1" -QUBOTools = "0.4" -julia = "1.6" +MathOptInterface = "1" +QUBOTools = "0.6.1" +julia = "1.6" diff --git a/docs/Project.toml b/docs/Project.toml index b6370474..96713b1e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,18 +1,20 @@ [deps] Anneal = "e4d9eb7f-b088-426e-aeb5-1c0dae3d8abb" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +DWaveNeal = "870cdf72-5502-4b10-839c-127ceab78f22" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +MQLib = "16f11440-1623-44c9-850c-358a6c72f3c9" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" ToQUBO = "9a412ddf-83fa-43b6-9748-7843c851aa65" [compat] -Anneal = "0.1" +Anneal = "0.6" CSV = "0.10" DataFrames = "1.3" Documenter = "0.27" -JuMP = "0.23" +JuMP = "1" MathOptInterface = "1" ToQUBO = "0.1" diff --git a/docs/make.jl b/docs/make.jl index 0ee002cf..c9832a7d 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -19,14 +19,21 @@ makedocs(; authors="Pedro Xavier and Tiago Andrade and Joaquim Garcia and David Bernal", pages=[ "Home" => "index.md", - "manual.md", - "examples.md", + "Manual" => "manual.md", + "Examples" => [ + "Knapsack" =>"examples/knapsack.md", + "Prime Factorization" => "examples/prime_factorization.md", + ], "Booklet" => "booklet.md" ], workdir="." ) -deploydocs( - repo=raw"github.com/psrenergy/ToQUBO.jl.git", - push_preview = true -) \ No newline at end of file +if "--skip-deploy" ∈ ARGS + @warn "Skipping deployment" +else + deploydocs( + repo=raw"github.com/psrenergy/ToQUBO.jl.git", + push_preview = true + ) +end \ No newline at end of file diff --git a/docs/src/booklet.md b/docs/src/booklet.md index 3c3161d1..5eec2071 100644 --- a/docs/src/booklet.md +++ b/docs/src/booklet.md @@ -18,8 +18,12 @@ ToQUBO.toqubo! ```@docs ToQUBO.toqubo_sense! ToQUBO.toqubo_variables! -ToQUBO.toqubo_constraint! +ToQUBO.toqubo_constraint +ToQUBO.toqubo_constraints! ToQUBO.toqubo_objective! +ToQUBO.toqubo_penalties! +ToQUBO.toqubo_parse! +ToQUBO.toqubo_build! ``` ## Pseudo-Boolean Optimization @@ -32,16 +36,23 @@ ToQUBO.PBO.derivative ToQUBO.PBO.gradient ToQUBO.PBO.gap ToQUBO.PBO.sharpness -ToQUBO.PBO.discretize ToQUBO.PBO.relaxed_gcd ``` -### Quadratization +### A Primer on Submodularity +A set function ``f : 2^{S} \to \mathbb{R}`` is said to be submodular if + +```math +f(X \cup Y) + f(X \cap Y) \le f(X) + f(Y) \forall X, Y \subset S +``` + +holds. + +## Quadratization In order to successfully achieve a QUBO formulation, sometimes it is needed to quadratize the resulting PBF, i.e., reduce its degree until reaching the quadratic case. There are many quadratization methods available, and `ToQUBO` implements a few of them. ```@docs -ToQUBO.PBO.quadratize -ToQUBO.PBO.@quadratization +ToQUBO.PBO.quadratize! ``` ## Virtual Mapping @@ -49,24 +60,39 @@ During reformulation, `ToQUBO` holds two distinct models, namely the *Source Mod This is done in a transparent fashion for both agents since the user will mostly interact with the presented model, and the solvers will only access the generated one. -### Virtual Variables +## Virtual Variables Every virtual model stores a collection of virtual variables, intended to provide a link between those in the source and those to be created in the target model. Each virtual variable stores enconding information for later expansion and evaluation. ```@docs -ToQUBO.VirtualMapping.VirtualVariable -ToQUBO.VirtualMapping.mapvar! -ToQUBO.VirtualMapping.expandℝ! -ToQUBO.VirtualMapping.slackℝ! -ToQUBO.VirtualMapping.expandℤ! -ToQUBO.VirtualMapping.slackℤ! -ToQUBO.VirtualMapping.mirror𝔹! -ToQUBO.VirtualMapping.slack𝔹! +ToQUBO.VirtualVariable +ToQUBO.encode! +``` + +## Variable Encoding + +### Linear Encoding Methods +```@docs +ToQUBO.LinearEncoding +ToQUBO.Binary +ToQUBO.Unary +ToQUBO.Arithmetic +ToQUBO.OneHot +``` + +### Sequential Encoding Methods +```@docs +ToQUBO.SequentialEncoding +ToQUBO.DomainWall +``` + +### Bounded Coefficients +```@docs +ToQUBO.Bounded ``` -### Virtual Models +## Virtual Models ```@docs -ToQUBO.VirtualMapping.AbstractVirtualModel -ToQUBO.VirtualQUBOModel +ToQUBO.VirtualModel ``` ### Annealing & Sampling diff --git a/docs/src/examples.md b/docs/src/examples/knapsack.md similarity index 78% rename from docs/src/examples.md rename to docs/src/examples/knapsack.md index 778abf56..3a40ecd1 100644 --- a/docs/src/examples.md +++ b/docs/src/examples/knapsack.md @@ -1,5 +1,3 @@ -# Examples - ## Knapsack We start with some instances of the discrete [Knapsack Problem](https://en.wikipedia.org/wiki/Knapsack_problem) whose standard formulation is ```math @@ -18,13 +16,13 @@ import MathOptInterface as MOI const MOIU = MOI.Utilities using ToQUBO -using Anneal # <- Your favourite Annealer / Sampler / Solver here +using DWaveNeal # <- Your favourite Annealer / Sampler / Solver here # Example from https://jump.dev/MathOptInterface.jl/stable/tutorials/example/ # Virtual QUBO Model model = MOI.instantiate( - () -> ToQUBO.Optimizer(SimulatedAnnealer.Optimizer), + () -> ToQUBO.Optimizer(DWaveNeal.Optimizer), with_bridge_type = Float64, ) @@ -33,10 +31,12 @@ c = [1.0, 2.0, 3.0] w = [0.3, 0.5, 1.0] C = 3.2; -# -*- Variables -*- # x = MOI.add_variables(model, n); -# -*- Objective -*- # +for xᵢ in x + MOI.add_constraint(model, xᵢ, MOI.ZeroOne()) +end + MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) MOI.set( @@ -45,22 +45,16 @@ MOI.set( MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, x), 0.0), ); -# -*- Constraints -*- # MOI.add_constraint( model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(w, x), 0.0), MOI.LessThan(C), ); -for xᵢ in x - MOI.add_constraint(model, xᵢ, MOI.ZeroOne()) -end - -# Run! MOI.optimize!(model) # Collect Solution -MOI.get(model, MOI.VariablePrimal(), x) +MOI.get.(model, MOI.VariablePrimal(), x) ``` ### JuMP + D-Wave Examples @@ -73,7 +67,7 @@ using CSV using DataFrames using Random -# -> Generate Data <- +# Generate Data rng = MersenneTwister(1) df = DataFrame( @@ -94,35 +88,21 @@ df = CSV.read("knapsack.csv", DataFrame) ```@example dwave-knapsack using JuMP using ToQUBO -using Anneal # <- Your favourite Annealer / Sampler / Solver here +using DWaveNeal # <- Your favourite Annealer/Sampler/Solver here -# -> Model <- -model = Model(() -> ToQUBO.Optimizer(SimulatedAnnealer.Optimizer)) +model = Model(() -> ToQUBO.Optimizer(DWaveNeal.Optimizer)) n = size(df, 1) c = collect(Float64, df[!, :cost]) w = collect(Float64, df[!, :weight]) C = round(0.8 * sum(w)) -# -> Variables <- -@variable(model, x[i=1:n], Bin) - -# -> Objective <- +@variable(model, x[1:n], Bin) @objective(model, Max, c' * x) - -# -> Constraint <- @constraint(model, w' * x <= C) -# ->-> Run! ->-> optimize!(model) # Add Results as a new column -insertcols!( - df, - 3, - :select => map( - (ξ) -> (ξ > 0.0) ? "Yes" : "No", - value.(x), - ), -) +insertcols!(df, 3, :select => map((ξ) -> (ξ > 0.0) ? "Yes" : "No", value.(x))) ``` \ No newline at end of file diff --git a/docs/src/examples/prime_factoring.md b/docs/src/examples/prime_factoring.md deleted file mode 100644 index 46161570..00000000 --- a/docs/src/examples/prime_factoring.md +++ /dev/null @@ -1,21 +0,0 @@ -# Prime Factoring - -```@setup prime-factoring -using JuMP -using ToQUBO -``` - -```@example prime-factoring -P = 3 -Q = 5 -R = P * Q - -model = Model(() -> ToQUBO.Optimizer(ExactSampler.Optimizer)) - -@variable(model, 0 <= p <= P, Integer) -@variable(model, 0 <= q <= Q, Integer) - -@objective(model, Min, (R - p * q) ^ 2) - -optimize!(model) -``` \ No newline at end of file diff --git a/docs/src/examples/prime_factorization.md b/docs/src/examples/prime_factorization.md new file mode 100644 index 00000000..86391dce --- /dev/null +++ b/docs/src/examples/prime_factorization.md @@ -0,0 +1,57 @@ +# Prime Factorization + +A central problem in Number Theory and cryptography is to factor ``R \in \mathbb{N}``, which is known +to be the product of two distinct prime numbers ``p, q \in \mathbb{N}``. +[Shor's Algorithm](https://en.wikipedia.org/wiki/Shor%27s_algorithm), designed to address such task is +often regarded as one of the major theoretical landmarks in Quantum Computing, being responsible for +driving increasingly greater interest to the area. + +A naïve approach to model this problem can be stated as a quadratically-constrained integer program: +```math +\begin{array}{rl} +\text{s.t.} & p \times q = R \\ + & p, q \ge 0 \\ + & p, q \in \mathbb{Z} +\end{array} +``` + +From the definition and the basics of number theory, we are able to retrieve a few assumptions about the problem's variables: +- ``p`` and ``q`` are bounded to the interval ``\left[1, R\right]`` +- Moreover, it is fine to assume ``1 < p \le \sqrt{R} \le q \le R \div 2``. + +```@example prime-factorization +using JuMP +using ToQUBO +using DWaveNeal + +function factor(R::Integer; optimizer = DWaveNeal.Optimizer) + return factor(identity, R; optimizer = optimizer) +end + +function factor(config!::Function, R::Integer; optimizer = DWaveNeal.Optimizer) + model = Model(() -> ToQUBO.Optimizer(optimizer)) + + @variable(model, 1 <= p <= √R, Int) + @variable(model, √R <= q <= R ÷ 2, Int) + + @constraint(model, p * q == R) + + config!(model) + + optimize!(model) + + p = round(Int, value(model[:p])) + q = round(Int, value(model[:q])) + + return (p, q) +end +``` + +```@example prime-factorization +p, q = factor(15) do model + set_optimizer_attribute(model, "num_reads", 1_000) + set_optimizer_attribute(model, "num_sweeps", 2_000) +end + +print("$p, $q") +``` \ No newline at end of file diff --git a/docs/src/manual.md b/docs/src/manual.md index c0ac5905..bc4c311a 100644 --- a/docs/src/manual.md +++ b/docs/src/manual.md @@ -4,15 +4,38 @@ ```@example quick-start using JuMP using ToQUBO -using Anneal +using DWaveNeal # <- Your favourite Annealer/Sampler/Solver here -model = Model(() -> ToQUBO.Optimizer(SimulatedAnnealer.Optimizer)) +model = Model(() -> ToQUBO.Optimizer(DWaveNeal.Optimizer)) @variable(model, x[1:3], Bin) + @objective(model, Max, 1.0 * x[1] + 2.0 * x[2] + 3.0 * x[3]) -@constraint(model, 0.3 * x[1] + 0.5 * x[2] + 1.0 * x[3] <= 3.2) + +@constraint(model, 0.3 * x[1] + 0.5 * x[2] + 1.0 * x[3] <= 1.6) optimize!(model) solution_summary(model) +``` + +## Compiler Flags + +### Architecture +```@docs +ToQUBO.Attributes.Architecture +``` + +### Quadratization +```@docs +ToQUBO.Attributes.Quadratize +ToQUBO.Attributes.QuadratizationMethod +ToQUBO.Attributes.StableQuadratization +``` + +### Variable & Constraint Encoding +```@docs +ToQUBO.Attributes.VariableEncodingMethod +ToQUBO.Attributes.VariableEncodingPenalty +ToQUBO.Attributes.ConstraintEncodingPenalty ``` \ No newline at end of file diff --git a/src/ToQUBO.jl b/src/ToQUBO.jl index 77f450c7..a0d6de35 100644 --- a/src/ToQUBO.jl +++ b/src/ToQUBO.jl @@ -1,18 +1,21 @@ module ToQUBO -# -*- :: Base Imports & Constants :: -*- # -using TOML -using Base: @kwdef +# Base Imports & Constants +import TOML const PROJECT_FILE_PATH = joinpath(@__DIR__, "..", "Project.toml") const PROJECT_VERSION = VersionNumber(getindex(TOML.parsefile(PROJECT_FILE_PATH), "version")) -# -*- :: External Imports :: -*- # -using MathOptInterface -const MOI = MathOptInterface +# QUBOTools +import QUBOTools + +const QUBO_NORMAL_FORM{T} = Tuple{Int,Dict{Int,T},Dict{Tuple{Int,Int},T},T,T} + +# External Imports +import MathOptInterface as MOI const MOIU = MOI.Utilities const MOIB = MOI.Bridges -# -*- MOI Aliases -*- # +# MOI Aliases const SAF{T} = MOI.ScalarAffineFunction{T} const SQF{T} = MOI.ScalarQuadraticFunction{T} const SAT{T} = MOI.ScalarAffineTerm{T} @@ -25,23 +28,23 @@ const GT{T} = MOI.GreaterThan{T} const VI = MOI.VariableIndex const CI = MOI.ConstraintIndex -# -*- :: QUBOTools :: -*- # -import QUBOTools: QUBOTools, qubo, backend - -# -*- :: Library Icludes :: -*- # - -# -*- Library -*- # +# Library include("lib/error.jl") include("lib/pbo/PBO.jl") -# -*- Model -*- # +# Model include("model/qubo.jl") include("model/prequbo.jl") include("model/virtual.jl") include("model/wrapper.jl") -include("model/attributes.jl") -# ~*~ Compiler & Analysis ~*~ # +# Compiler & Analysis include("compiler/compiler.jl") +# Attributes +include("attributes/model.jl") +include("attributes/solver.jl") +include("attributes/virtual.jl") +include("attributes/compiler.jl") + end # module \ No newline at end of file diff --git a/src/attributes/compiler.jl b/src/attributes/compiler.jl new file mode 100644 index 00000000..b0396078 --- /dev/null +++ b/src/attributes/compiler.jl @@ -0,0 +1,341 @@ +module Attributes + +import ..ToQUBO +import MathOptInterface as MOI +const MOIU = MOI.Utilities +const VI = MOI.VariableIndex +const CI = MOI.ConstraintIndex + +abstract type CompilerAttribute <: MOI.AbstractOptimizerAttribute end + +export + Architecture, + Discretize, + Quadratize, + QuadratizationMethod, + StableQuadratization, + DefaultVariableEncodingATol, + DefaultVariableEncodingBits, + DefaultVariableEncodingMethod, + VariableEncodingATol, + VariableEncodingBits, + VariableEncodingMethod, + VariableEncodingPenalty, + ConstraintEncodingPenalty, + QUBONormalForm + +@doc raw""" + Architecture() +""" struct Architecture <: CompilerAttribute end + +function MOI.get(model::ToQUBO.VirtualModel, ::Architecture)::ToQUBO.AbstractArchitecture + return get(model.compiler_settings, :architecture, ToQUBO.GenericArchitecture()) +end + +function MOI.set(model::ToQUBO.VirtualModel, ::Architecture, arch::ToQUBO.AbstractArchitecture) + model.compiler_settings[:architecture] = arch + + return nothing +end + +@doc raw""" + Discretize() + +When set, this boolean flag guarantees that every coefficient in the final formulation is an integer. +""" struct Discretize <: CompilerAttribute end + +function MOI.get(model::ToQUBO.VirtualModel, ::Discretize, flag::Bool)::Bool + return get(model.compiler_settings, :discretize, false) +end + +function MOI.set(model::ToQUBO.VirtualModel, ::Discretize, flag::Bool) + model.compiler_settings[:discretize] = flag + + return nothing +end + +@doc raw""" + Quadratize() + +Boolean flag to conditionally perform the quadratization step. +Is automatically set by the compiler when high-order functions are generated. +""" struct Quadratize <: CompilerAttribute end + +function MOI.get(model::ToQUBO.VirtualModel, ::Quadratize)::Bool + return get(model.compiler_settings, :quadratize, false) +end + +function MOI.set(model::ToQUBO.VirtualModel, ::Quadratize, flag::Bool) + model.compiler_settings[:quadratize] = flag + + return nothing +end + +@doc raw""" + QuadratizationMethod() + +Defines which quadratization method to use. +Available options are defined in the [`PBO`](@ref) submodule. +""" struct QuadratizationMethod <: CompilerAttribute end + +function MOI.get(model::ToQUBO.VirtualModel, ::QuadratizationMethod) + return get(model.compiler_settings, :QuadratizationMethod, ToQUBO.PBO.INFER) +end + +function MOI.set( + model::ToQUBO.VirtualModel, + ::QuadratizationMethod, + ::Type{method}, +) where {method<:ToQUBO.PBO.QuadratizationMethod} + model.compiler_settings[:QuadratizationMethod] = method + + return nothing +end + +@doc raw""" + StableQuadratization() + +When set, this boolean flag enables stable quadratization methods, thus yielding predictable results. +This is intended to be used during tests or other situations where deterministic output is desired. +On the other hand, usage in production is not recommended since it requires increased memory and processing resources. +""" struct StableQuadratization <: CompilerAttribute end + +function MOI.get(model::ToQUBO.VirtualModel, ::StableQuadratization)::Bool + return get(model.compiler_settings, :stable_quadratization, false) +end + +function MOI.set(model::ToQUBO.VirtualModel, ::StableQuadratization, flag::Bool) + model.compiler_settings[:stable_quadratization] = flag + + return nothing +end + +@doc raw""" + DefaultVariableEncodingMethod() +""" struct DefaultVariableEncodingMethod <: CompilerAttribute end + +function MOI.get(model::ToQUBO.VirtualModel, ::DefaultVariableEncodingMethod)::ToQUBO.Encoding + return get(model.compiler_settings, :default_variable_encoding_method, ToQUBO.Binary()) +end + +function MOI.set(model::ToQUBO.VirtualModel, ::DefaultVariableEncodingMethod, e::ToQUBO.Encoding) + model.compiler_settings[:default_variable_encoding_method] = e + + return nothing +end + +@doc raw""" + DefaultVariableEncodingATol() +""" struct DefaultVariableEncodingATol <: CompilerAttribute end + +function MOI.get(model::ToQUBO.VirtualModel{T}, ::DefaultVariableEncodingATol)::T where {T} + return get(model.compiler_settings, :default_variable_encoding_atol, T(1E-2)) +end + +function MOI.set(model::ToQUBO.VirtualModel{T}, ::DefaultVariableEncodingATol, τ::T) where {T} + model.compiler_settings[:default_variable_encoding_atol] = τ + + return nothing +end + +abstract type CompilerVariableAttribute <: MOI.AbstractVariableAttribute end + +@doc raw""" + VariableEncodingATol() +""" struct VariableEncodingATol <: CompilerVariableAttribute end + +function MOI.get(model::ToQUBO.VirtualModel{T}, ::VariableEncodingATol, vi::VI)::T where {T} + attr = :variable_encoding_atol + + if !haskey(model.variable_settings, attr) || !haskey(model.variable_settings[attr], vi) + return MOI.get(model, DefaultVariableEncodingATol()) + else + return model.variable_settings[attr][vi] + end +end + +function MOI.set(model::ToQUBO.VirtualModel{T}, ::VariableEncodingATol, vi::VI, τ::T) where {T} + attr = :variable_encoding_atol + + if !haskey(model.variable_settings, attr) + model.variable_settings[attr] = Dict{VI,Any}(vi => τ) + else + model.variable_settings[attr][vi] = τ + end + + return nothing +end + +@doc raw""" + DefaultVariableEncodingBits() +""" struct DefaultVariableEncodingBits <: CompilerAttribute end + +function MOI.get(model::ToQUBO.VirtualModel, ::DefaultVariableEncodingBits)::Union{Integer,Nothing} + return get(model.compiler_settings, :default_variable_encoding_bits, nothing) +end + +function MOI.set(model::ToQUBO.VirtualModel, ::DefaultVariableEncodingBits, n::Union{Integer,Nothing}) + model.compiler_settings[:default_variable_encoding_bits] = n + + return nothing +end + +@doc raw""" + VariableEncodingBits() +""" struct VariableEncodingBits <: CompilerVariableAttribute end + +function MOI.get(model::ToQUBO.VirtualModel, ::VariableEncodingBits, vi::VI)::Union{Integer,Nothing} + attr = :variable_encoding_bits + + if !haskey(model.variable_settings, attr) || !haskey(model.variable_settings[attr], vi) + return MOI.get(model, DefaultVariableEncodingBits()) + else + return model.variable_settings[attr][vi] + end +end + +function MOI.set(model::ToQUBO.VirtualModel, ::VariableEncodingBits, vi::VI, n::Union{Integer,Nothing}) + attr = :variable_encoding_bits + + if !haskey(model.variable_settings, attr) + model.variable_settings[attr] = Dict{VI,Any}(vi => n) + else + model.variable_settings[attr][vi] = n + end + + return nothing +end + +@doc raw""" + VariableEncodingMethod() + +Available methods are: +- [`Binary`](@ref) (default) +- [`Unary`](@ref) +- [`Arithmetic`](@ref) +- [`OneHot`](@ref) +- [`DomainWall`](@ref) +- [`Bounded`](@ref) + +The [`Binary`](@ref), [`Unary`](@ref) and [`Arithmetic`](@ref) encodings can have their +expansion coefficients bounded by parametrizing the [`Bounded`](@ref) encoding. +""" struct VariableEncodingMethod <: CompilerVariableAttribute end + +function MOI.get(model::ToQUBO.VirtualModel, ::VariableEncodingMethod, vi::VI)::ToQUBO.Encoding + attr = :variable_encoding_method + + if !haskey(model.variable_settings, attr) || !haskey(model.variable_settings[attr], vi) + return MOI.get(model, DefaultVariableEncodingMethod()) + else + return model.variable_settings[attr][vi] + end +end + +function MOI.set(model::ToQUBO.VirtualModel, ::VariableEncodingMethod, vi::VI, e::ToQUBO.Encoding) + attr = :variable_encoding_method + + if !haskey(model.variable_settings, attr) + model.variable_settings[attr] = Dict{VI,Any}(vi => e) + else + model.variable_settings[attr][vi] = e + end + + return nothing +end + +@doc raw""" + VariableEncodingPenalty() + +Allows the user to set and retrieve the coefficients used for encoding variables when additional +constraints are involved. +""" struct VariableEncodingPenalty <: CompilerVariableAttribute end + +function MOI.get(model::ToQUBO.VirtualModel{T}, ::VariableEncodingPenalty, vi::VI) where {T} + return model.θ[vi]::T +end + +function MOI.set( + model::ToQUBO.VirtualModel{T}, + ::VariableEncodingPenalty, + vi::VI, + θ::T, +) where {T} + model.θ[vi] = θ + + return nothing +end + +MOI.is_set_by_optimize(::VariableEncodingPenalty) = true + +abstract type CompilerConstraintAttribute <: MOI.AbstractConstraintAttribute end + +@doc raw""" + ConstraintEncodingPenalty() + +Allows the user to set and retrieve the coefficients used for encoding constraints. +""" struct ConstraintEncodingPenalty <: CompilerConstraintAttribute end + +function MOI.get(model::ToQUBO.VirtualModel{T}, ::ConstraintEncodingPenalty, ci::CI) where {T} + return model.ρ[ci] +end + +function MOI.set( + model::ToQUBO.VirtualModel{T}, + ::ConstraintEncodingPenalty, + ci::CI, + ρ::T, +) where {T} + model.ρ[ci] = ρ + + return nothing +end + +MOI.is_set_by_optimize(::ConstraintEncodingPenalty) = true + +@doc raw""" + QUBONormalForm +""" struct QUBONormalForm <: CompilerAttribute end + +function MOI.get( + model::ToQUBO.VirtualModel{T}, + ::QUBONormalForm, +)::ToQUBO.QUBO_NORMAL_FORM{T} where {T} + target_model = model.target_model + + n = MOI.get(target_model, MOI.NumberOfVariables()) + F = MOI.get(target_model, MOI.ObjectiveFunctionType()) + f = MOI.get(target_model, MOI.ObjectiveFunction{F}()) + + linear_terms = sizehint!(Dict{Int,T}(), length(f.affine_terms)) + quadratic_terms = sizehint!(Dict{Tuple{Int,Int},T}(), length(f.quadratic_terms)) + + for a in f.affine_terms + c = a.coefficient + i = a.variable.value + + linear_terms[i] = get(linear_terms, i, zero(T)) + c + end + + for q in f.quadratic_terms + c = q.coefficient + i = q.variable_1.value + j = q.variable_2.value + + if i == j + linear_terms[i] = get(linear_terms, i, zero(T)) + c / 2 + elseif i > j + quadratic_terms[(j, i)] = get(quadratic_terms, (j, i), zero(T)) + c + else + quadratic_terms[(i, j)] = get(quadratic_terms, (i, j), zero(T)) + c + end + end + + scale = one(T) + offset = f.constant + + return (n, linear_terms, quadratic_terms, scale, offset) +end + +MOIU.map_indices(::MOIU.IndexMap, x::ToQUBO.QUBO_NORMAL_FORM{T}) where {T} = x + +end # module Settings \ No newline at end of file diff --git a/src/attributes/model.jl b/src/attributes/model.jl new file mode 100644 index 00000000..044309df --- /dev/null +++ b/src/attributes/model.jl @@ -0,0 +1,74 @@ +MOI.get(::VirtualModel, ::MOI.SolverName) = "Virtual Model" +MOI.get(::VirtualModel, ::MOI.SolverVersion) = PROJECT_VERSION + +const MOI_MODEL_ATTRIBUTE = Union{ + MOI.ListOfConstraintAttributesSet, + MOI.ListOfConstraintIndices, + MOI.ListOfConstraintTypesPresent, + MOI.ListOfModelAttributesSet, + MOI.ListOfVariableAttributesSet, + MOI.ListOfVariableIndices, + MOI.NumberOfConstraints, + MOI.NumberOfVariables, + MOI.Name, + MOI.ObjectiveFunction, + MOI.ObjectiveFunctionType, + MOI.ObjectiveSense, +} + +function MOI.get(model::VirtualModel, attr::MOI_MODEL_ATTRIBUTE) + return MOI.get(model.source_model, attr) +end + +function MOI.set(model::VirtualModel, attr::MOI_MODEL_ATTRIBUTE, value::Any) + MOI.set(model.source_model, attr, value) + + return nothing +end + +function MOI.get( + model::VirtualModel, + attr::Union{MOI.ConstraintFunction,MOI.ConstraintSet}, + ci::MOI.ConstraintIndex, +) + return MOI.get(model.source_model, attr, ci) +end + +function MOI.get(model::VirtualModel, attr::MOI.VariableName, x::VI) + return MOI.get(model.source_model, attr, x) +end + +function Base.show(io::IO, model::VirtualModel) + print( + io, + """ + Virtual Model + with source: + $(model.source_model) + with target: + $(model.target_model) + """, + ) +end + +function MOI.add_variable(model::VirtualModel) + return MOI.add_variable(model.source_model) +end + +function MOI.add_constraint( + model::VirtualModel, + f::MOI.AbstractFunction, + s::MOI.AbstractSet, +) + return MOI.add_constraint(model.source_model, f, s) +end + +function MOI.set( + model::VirtualModel, + ::MOI.ObjectiveFunction{F}, + f::F, +) where {F<:MOI.AbstractFunction} + MOI.set(model.source_model, MOI.ObjectiveFunction{F}(), f) + + return nothing +end \ No newline at end of file diff --git a/src/attributes/solver.jl b/src/attributes/solver.jl new file mode 100644 index 00000000..fbb2c2b1 --- /dev/null +++ b/src/attributes/solver.jl @@ -0,0 +1,110 @@ +function MOI.get(model::VirtualModel, attr::MOI.AbstractOptimizerAttribute) + if !isnothing(model.optimizer) && MOI.supports(model.optimizer, attr) + return MOI.get(model.optimizer, attr) + else + return MOI.get(model.source_model, attr) + end +end + +function MOI.set(model::VirtualModel, attr::MOI.AbstractOptimizerAttribute, value::Any) + if !isnothing(model.optimizer) && MOI.supports(model.optimizer, attr) + MOI.set(model.optimizer, attr, value) + else + MOI.set(model.source_model, attr, value) + end + + return nothing +end + +function MOI.supports(model::VirtualModel, attr::MOI.AbstractOptimizerAttribute) + if !isnothing(model.optimizer) + return MOI.supports(model.optimizer, attr) + else + return MOI.supports(model.source_model, attr) + end +end + +function MOI.get( + model::VirtualModel, + attr::Union{ + MOI.SolveTimeSec, + MOI.PrimalStatus, + MOI.DualStatus, + MOI.TerminationStatus, + MOI.RawStatusString, + }, +) + if !isnothing(model.optimizer) + return MOI.get(model.optimizer, attr) + else + return nothing + end +end + +function MOI.supports( + model::VirtualModel, + attr::Union{ + MOI.SolveTimeSec, + MOI.PrimalStatus, + MOI.DualStatus, + MOI.TerminationStatus, + MOI.RawStatusString, + }, +) + if !isnothing(model.optimizer) + return MOI.supports(model.optimizer, attr) + else + return false + end +end + +function MOI.get(model::VirtualModel, rc::MOI.ResultCount) + if isnothing(model.optimizer) + return 0 + else + return MOI.get(model.optimizer, rc) + end +end + +MOI.supports(::VirtualModel, ::MOI.ResultCount) = true + +function MOI.get(model::VirtualModel{T}, ov::MOI.ObjectiveValue) where {T} + if isnothing(model.optimizer) + return zero(T) + else + return MOI.get(model.optimizer, ov) + end +end + +function MOI.get(model::VirtualModel{T}, vp::MOI.VariablePrimalStart, x::VI) where {T} + return MOI.get(model.source_model, vp, x) +end + +MOI.supports(::VirtualModel, ::MOI.VariablePrimalStart, ::MOI.VariableIndex) = true + +function MOI.get(model::VirtualModel{T}, vp::MOI.VariablePrimal, x::VI) where {T} + if isnothing(model.optimizer) + return zero(T) + else + s = zero(T) + v = model.source[x] + + for (ω, c) in expansion(v) + for y in ω + c *= MOI.get(model.optimizer, vp, y) + end + + s += c + end + + return s + end +end + +function MOI.get(model::VirtualModel, rs::MOI.RawSolver) + if isnothing(model.optimizer) + return nothing + else + return MOI.get(model.optimizer, rs) + end +end diff --git a/src/attributes/virtual.jl b/src/attributes/virtual.jl new file mode 100644 index 00000000..50ad497a --- /dev/null +++ b/src/attributes/virtual.jl @@ -0,0 +1,19 @@ +function MOI.is_empty(model::VirtualModel) + return MOI.is_empty(model.source_model) +end + +function MOI.empty!(model::VirtualModel) + # Source Model + MOI.empty!(model.source_model) + + # Underlying Optimizer + if !isnothing(model.optimizer) + MOI.empty!(model.optimizer) + end + + return nothing +end + +function MOI.get(::VirtualModel{T}, ::MOIB.ListOfNonstandardBridges{T}) where {T} + return Type[] +end diff --git a/src/compiler/analysis.jl b/src/compiler/analysis.jl new file mode 100644 index 00000000..4e614634 --- /dev/null +++ b/src/compiler/analysis.jl @@ -0,0 +1,50 @@ +function is_qubo(model::MOI.ModelLike) + return is_quadratic(model) && + is_unconstrained(model) && + is_binary(model) && + is_optimization(model) +end + +function is_quadratic(model::MOI.ModelLike) + return MOI.get(model, MOI.ObjectiveFunctionType()) <: Union{SQF,SAF,VI} +end + +function is_unconstrained(model::MOI.ModelLike) + for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent()) + if !(F === VI && S === MOI.ZeroOne) + return false + end + end + + return true +end + +function is_binary(model::MOI.ModelLike) + m = MOI.get(model, MOI.NumberOfConstraints{VI,MOI.ZeroOne}()) + n = MOI.get(model, MOI.NumberOfVariables()) + + # There will be uncovered variables in this case: + if m < n + return false + end + + 𝔹 = sizehint!(Set{VI}(), n) + + for ci in MOI.get(model, MOI.ListOfConstraintIndices{VI,MOI.ZeroOne}()) + push!(𝔹, MOI.get(model, MOI.ConstraintFunction(), ci)) + end + + for vi in MOI.get(model, MOI.ListOfVariableIndices()) + if vi ∉ 𝔹 # Non-binary variable found + return false + end + end + + return true +end + +function is_optimization(model::MOI.ModelLike) + S = MOI.get(model, MOI.ObjectiveSense()) + + return (S === MOI.MAX_SENSE || S === MOI.MIN_SENSE) +end diff --git a/src/compiler/architectures.jl b/src/compiler/architectures.jl index d15c6822..52c0cefe 100644 --- a/src/compiler/architectures.jl +++ b/src/compiler/architectures.jl @@ -1,8 +1,10 @@ -abstract type AbstractArchitecture end +@doc raw""" + AbstractArchitecture +""" abstract type AbstractArchitecture end struct GenericArchitecture <: AbstractArchitecture end @doc raw""" """ function infer_architecture end -infer_architecture(::Any) = GenericArchitecture() \ No newline at end of file +infer_architecture(::Any) = GenericArchitecture() diff --git a/src/compiler/build.jl b/src/compiler/build.jl index 13f8ce52..77d88723 100644 --- a/src/compiler/build.jl +++ b/src/compiler/build.jl @@ -1,8 +1,27 @@ -using MutableArithmetics -const MA = MutableArithmetics +function toqubo_build!(model::VirtualModel{T}, arch::AbstractArchitecture) where {T} + # Assemble Objective Function + toqubo_hamiltonian!(model, arch) + + # Quadratization Step + toqubo_quadratize!(model, arch) -function toqubo_hamiltonian!(model::VirtualQUBOModel{T}, ::AbstractArchitecture) where {T} - copy!(model.H, model.f) + # Write to MathOptInterface + toqubo_output!(model, arch) + + return nothing +end + +function toqubo_hamiltonian!(model::VirtualModel{T}, ::AbstractArchitecture) where {T} + empty!(model.H) + + # Calculate an upper bound on the number of terms + num_terms = length(model.f) + sum(length, model.g; init=0) + sum(length, model.h; init=0) + + sizehint!(model.H, num_terms) + + for (ω, c) in model.f + model.H[ω] += c + end for (ci, g) in model.g ρ = model.ρ[ci] @@ -23,18 +42,18 @@ function toqubo_hamiltonian!(model::VirtualQUBOModel{T}, ::AbstractArchitecture) return nothing end -function toqubo_aux(model::VirtualQUBOModel, ::Nothing, ::AbstractArchitecture) - target_model = MOI.get(model, TargetModel()) +function toqubo_aux(model::VirtualModel, ::Nothing, ::AbstractArchitecture)::VI + target_model = model.target_model w = MOI.add_variable(target_model) MOI.add_constraint(target_model, w, MOI.ZeroOne()) - return w::VI + return w end -function toqubo_aux(model::VirtualQUBOModel, n::Integer, ::AbstractArchitecture)::Vector{VI} - target_model = MOI.get(model, TargetModel()) +function toqubo_aux(model::VirtualModel, n::Integer, ::AbstractArchitecture)::Vector{VI} + target_model = model.target_model w = MOI.add_variables(target_model, n) @@ -43,17 +62,23 @@ function toqubo_aux(model::VirtualQUBOModel, n::Integer, ::AbstractArchitecture) return w end -function toqubo_quadratize!(model::VirtualQUBOModel, arch::AbstractArchitecture) - H = PBO.quadratize(model.H) do (n::Union{Integer,Nothing} = nothing) - return toqubo_aux(model, n, arch) - end +function toqubo_quadratize!(model::VirtualModel, arch::AbstractArchitecture) + if MOI.get(model, Attributes.Quadratize()) + method = MOI.get(model, Attributes.QuadratizationMethod()) + stable = MOI.get(model, Attributes.StableQuadratization()) - copy!(model.H, H) + PBO.quadratize!( + model.H, + PBO.Quadratization{method}(stable), + ) do (n::Union{Integer,Nothing} = nothing) + return toqubo_aux(model, n, arch) + end + end return nothing end -function toqubo_output!(model::VirtualQUBOModel{T}, ::AbstractArchitecture) where {T} +function toqubo_output!(model::VirtualModel{T}, ::AbstractArchitecture) where {T} Q = SQT{T}[] a = SAT{T}[] b = zero(T) @@ -62,36 +87,26 @@ function toqubo_output!(model::VirtualQUBOModel{T}, ::AbstractArchitecture) wher if isempty(ω) b += c elseif length(ω) == 1 - push!(a, SAT{T}(c, ω...)) + x, = ω + + push!(a, SAT{T}(c, x)) elseif length(ω) == 2 - push!(Q, SQT{T}(c, ω...)) + x, y = ω + + push!(Q, SQT{T}(c, x, y)) else # NOTE: This should never happen in production. # During implementation of new quadratization and constraint reformulation methods # higher degree terms might be introduced by mistake. That's why it's important to # have this condition here. + # HINT: When debugging this, a good place to start is to check if the 'Quadratize' + # flag is set or not. If missing, it should mean that some constraint might induce + # PBFs of higher degree without calling 'MOI.set(model, Quadratize(), true)'. throw(QUBOError("Quadratization failed")) end end - MOI.set( - MOI.get(model, TargetModel()), - MOI.ObjectiveFunction{SQF{T}}(), - SQF{T}(Q, a, b), - ) + MOI.set(model.target_model, MOI.ObjectiveFunction{SQF{T}}(), SQF{T}(Q, a, b)) return nothing end - -function toqubo_build!(model::VirtualQUBOModel{T}, arch::AbstractArchitecture) where {T} - # -*- Assemble Objective Function -*- - toqubo_hamiltonian!(model, arch) - - # -*- Quadratization Step -*- - toqubo_quadratize!(model, arch) - - # -*- Write to MathOptInterface -*- - toqubo_output!(model, arch) - - return nothing -end \ No newline at end of file diff --git a/src/compiler/compiler.jl b/src/compiler/compiler.jl index ef9d2e3a..3a9836fe 100644 --- a/src/compiler/compiler.jl +++ b/src/compiler/compiler.jl @@ -1,6 +1,6 @@ +include("analysis.jl") include("architectures.jl") include("interface.jl") -include("validation.jl") include("parse.jl") include("variables.jl") include("objective.jl") @@ -8,13 +8,12 @@ include("constraints.jl") include("penalties.jl") include("build.jl") -# -*- toqubo: MOI.ModelLike -> QUBO.Model -*- +# toqubo: MOI.ModelLike -> QUBO.Model toqubo( source::MOI.ModelLike, arch::Union{AbstractArchitecture,Nothing} = nothing, optimizer = nothing, -) = toqubo(Float64, source, arch; optimizer = optimizer) - +) = toqubo(Float64, source, arch; optimizer) function toqubo( ::Type{T}, @@ -22,7 +21,7 @@ function toqubo( arch::Union{AbstractArchitecture,Nothing} = nothing; optimizer = nothing, ) where {T} - model = VirtualQUBOModel{T}(optimizer) + model = VirtualModel{T}(optimizer) MOI.copy_to(model, source) @@ -36,29 +35,95 @@ function toqubo( end function toqubo!( - model::VirtualQUBOModel{T}, + model::VirtualModel{T}, arch::AbstractArchitecture = GenericArchitecture(), ) where {T} - # :: Objective Sense :: # + # Cleanup + toqubo_empty!(model, arch) + + if is_qubo(model.source_model) + toqubo_copy!(model, arch) + + return nothing + end + + toqubo_compile!(model, arch) + + return nothing +end + +function toqubo_copy!( + model::VirtualModel{T}, + ::AbstractArchitecture, +) where {T} + source_model = model.source_model + target_model = model.target_model + + # Map Variables + for vi in MOI.get(source_model, MOI.ListOfVariableIndices()) + encode!(model, Mirror(), vi) + end + + # Copy Objective Sense + s = MOI.get(source_model, MOI.ObjectiveSense()) + + MOI.set(target_model, MOI.ObjectiveSense(), s) + + # Copy Objective Function + F = MOI.get(source_model, MOI.ObjectiveFunctionType()) + f = MOI.get(source_model, MOI.ObjectiveFunction{F}()) + + MOI.set(target_model, MOI.ObjectiveFunction{F}(), f) + + return nothing +end + +function toqubo_compile!( + model::VirtualModel{T}, + arch::AbstractArchitecture = GenericArchitecture(), +) where {T} + # Objective Sense toqubo_sense!(model, arch) - # :: Problem Variables :: # + # Problem Variables toqubo_variables!(model, arch) - # :: Objective Analysis :: # + # Objective Analysis toqubo_objective!(model, arch) - # :: Add Regular Constraints :: # + # Add Regular Constraints toqubo_constraints!(model, arch) - # :: Add Encoding Constraints :: # + # Add Encoding Constraints toqubo_encoding_constraints!(model, arch) - # :: Compute penalties :: # + # Compute penalties toqubo_penalties!(model, arch) - # :: Build Final Model :: # + # Build Final Model toqubo_build!(model, arch) return nothing -end \ No newline at end of file +end + +function toqubo_empty!( + model::VirtualModel, + ::AbstractArchitecture = GenericArchitecture(), +) + # Model + MOI.empty!(model.target_model) + + # Virtual Variables + empty!(model.variables) + empty!(model.source) + empty!(model.target) + + # PBF/IR + empty!(model.f) + empty!(model.g) + empty!(model.h) + empty!(model.ρ) + empty!(model.θ) + + return nothing +end diff --git a/src/compiler/constraints.jl b/src/compiler/constraints.jl index 59fdec69..e2c435ac 100644 --- a/src/compiler/constraints.jl +++ b/src/compiler/constraints.jl @@ -1,4 +1,4 @@ -function toqubo_constraints!(model::VirtualQUBOModel, arch::AbstractArchitecture) +function toqubo_constraints!(model::VirtualModel, arch::AbstractArchitecture) for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent()) for ci in MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) f = MOI.get(model, MOI.ConstraintFunction(), ci) @@ -14,25 +14,60 @@ function toqubo_constraints!(model::VirtualQUBOModel, arch::AbstractArchitecture return nothing end +@doc raw""" + toqubo_constraint( + ::VirtualModel{T}, + ::VI, + ::Union{ + MOI.ZeroOne, + MOI.Integer, + MOI.Interval{T}, + MOI.LessThan{T}, + MOI.GreaterThan{T} + }, + ::AbstractArchitecture + ) where {T} + +This method skips bound constraints over variables. +""" function toqubo_constraint( - ::VirtualQUBOModel{T}, + ::VirtualModel{T}, ::VI, ::Union{MOI.ZeroOne,MOI.Integer,MOI.Interval{T},LT{T},GT{T}}, ::AbstractArchitecture, -) where {T} end +) where {T} + return nothing +end + +@doc raw""" + toqubo_constraint(model::VirtualModel{T}, f::SAF{T}, s::EQ{T}, ::AbstractArchitecture) where {T} + +Turns constraints of the form + +```math +\begin{array}{rl} +\text{s.t} & \mathbf{a}'\mathbf{x} - b = 0 +\end{array} +``` + +into +```math +\left\Vertg(\mathbf{x})\right\Vert_{\left\lbrace{0}\right\rbrace} = \left(\mathbf{a}'\mathbf{x} - b\right)^{2} +``` +""" function toqubo_constraint( - model::VirtualQUBOModel{T}, + model::VirtualModel{T}, f::SAF{T}, s::EQ{T}, arch::AbstractArchitecture, ) where {T} - # -*- Scalar Affine Function: g(x) = a'x - b = 0 ~ 😄 -*- + # Scalar Affine Equality Constraint: g(x) = a'x - b = 0 g = toqubo_parse(model, f, s, arch) - + PBO.discretize!(g) - # -*- Bounds & Slack Variable -*- + # Bounds & Slack Variable l, u = PBO.bounds(g) if u < zero(T) # Always feasible @@ -45,18 +80,37 @@ function toqubo_constraint( return g^2 end +@doc raw""" + toqubo_constraint(model::VirtualModel{T}, f::SAF{T}, s::LT{T}, ::AbstractArchitecture) where {T} + +Turns constraints of the form + +```math +\begin{array}{rl} +\text{s.t} & \mathbf{a}'\mathbf{x} - b \le 0 +\end{array} +``` + +into + +```math +\left\Vertg(\mathbf{x})\right\Vert_{\left\lbrace{0}\right\rbrace} = \left(\mathbf{a}'\mathbf{x} - b\right + z)^{2} +``` + +by adding a slack variable ``z``. +""" function toqubo_constraint( - model::VirtualQUBOModel{T}, + model::VirtualModel{T}, f::SAF{T}, s::LT{T}, arch::AbstractArchitecture, ) where {T} - # -*- Scalar Affine Function: g(x) = a'x - b ≤ 0 🤔 -*- + # Scalar Affine Inequality Constraint: g(x) = a'x - b ≤ 0 g = toqubo_parse(model, f, s, arch) PBO.discretize!(g) - # -*- Bounds & Slack Variable -*- + # Bounds & Slack Variable l, u = PBO.bounds(g) if u < zero(T) # Always feasible @@ -67,7 +121,8 @@ function toqubo_constraint( end # Slack Variable - z = encode!(Binary(), model, nothing, zero(T), abs(l)) + e = MOI.get(model, Attributes.DefaultVariableEncodingMethod()) + z = encode!(model, e, nothing, zero(T), abs(l)) for (ω, c) in expansion(z) g[ω] += c @@ -77,17 +132,17 @@ function toqubo_constraint( end function toqubo_constraint( - model::VirtualQUBOModel{T}, + model::VirtualModel{T}, f::SQF{T}, s::EQ{T}, arch::AbstractArchitecture, ) where {T} - # -*- Scalar Quadratic Function: g(x) = x Q x + a x - b = 0 😢 -*- + # Scalar Quadratic Equality Constraint: g(x) = x Q x + a x - b = 0 g = toqubo_parse(model, f, s, arch) PBO.discretize!(g) - # -*- Bounds & Slack Variable -*- + # Bounds & Slack Variable l, u = PBO.bounds(g) if u < zero(T) # Always feasible @@ -97,21 +152,24 @@ function toqubo_constraint( @warn "Infeasible constraint detected" end + # Tell the compiler that quadratization is necessary + MOI.set(model, Attributes.Quadratize(), true) + return g^2 end function toqubo_constraint( - model::VirtualQUBOModel{T}, + model::VirtualModel{T}, f::SQF{T}, s::LT{T}, arch::AbstractArchitecture, ) where {T} - # -*- Scalar Quadratic Function: g(x) = x Q x + a x - b ≤ 0 😢 -*- + # Scalar Quadratic Inequality Constraint: g(x) = x Q x + a x - b ≤ 0 g = toqubo_parse(model, f, s, arch) PBO.discretize!(g) - # -*- Bounds & Slack Variable -*- + # Bounds & Slack Variable l, u = PBO.bounds(g) if u < zero(T) # Always feasible @@ -121,40 +179,57 @@ function toqubo_constraint( @warn "Infeasible constraint detected" end - z = encode!(Binary(), model, nothing, zero(T), abs(l)) + # Slack Variable + e = MOI.get(model, Attributes.DefaultVariableEncodingMethod()) + z = encode!(model, e, nothing, zero(T), abs(l)) for (ω, c) in expansion(z) g[ω] += c end + # Tell the compiler that quadratization is necessary + MOI.set(model, Attributes.Quadratize(), true) + return g^2 end function toqubo_constraint( - model::VirtualQUBOModel{T}, - v::MOI.VectorOfVariables, + model::VirtualModel{T}, + x::MOI.VectorOfVariables, ::MOI.SOS1{T}, ::AbstractArchitecture, ) where {T} - # -*- Special Ordered Set of Type 1: ∑ x <= min x 😄 -*- + # Special Ordered Set of Type 1: ∑ x ≤ min x g = PBO.PBF{VI,T}() - for vi in v.variables - for (ωi, _) in expansion(MOI.get(model, Source(), vi)) + for xi in x.variables + vi = model.source[xi] + + @assert encoding(vi) isa Mirror "Currently, SOS1 only supports binary variables" + + for (ωi, _) in expansion(vi) g[ωi] = one(T) end end - z = expansion(encode!(Mirror(), model, nothing)) + # Slack variable + e = Mirror() + z = encode!(model, e, nothing) - return (g + z - one(T))^2 # one-hot approach + for (ω, c) in expansion(z) + g[ω] += c + end + + g[nothing] += -one(T) + + return g^2 end function toqubo_encoding_constraints!( - model::VirtualQUBOModel{T}, + model::VirtualModel{T}, ::AbstractArchitecture, ) where {T} - for v in MOI.get(model, Variables()) + for v in model.variables if is_aux(v) continue end diff --git a/src/compiler/interface.jl b/src/compiler/interface.jl index 24a6a373..6a1d09df 100644 --- a/src/compiler/interface.jl +++ b/src/compiler/interface.jl @@ -18,74 +18,70 @@ For it to be true, a few conditions must be met: optimizer::Union{Nothing, Type{<:MOI.AbstractOptimizer}} = nothing ) -Low-level interface to create a `::VirtualQUBOModel{T}` from `::MOI.ModelLike` instance. +Low-level interface to create a `::VirtualModel{T}` from `::MOI.ModelLike` instance. If provided, an `::MOI.AbstractOptimizer` is attached to the model. """ function toqubo end @doc raw""" - toqubo!(model::VirtualQUBOModel{T}, ::AbstractArchitecture) where {T} + toqubo!(model::VirtualModel{T}, ::AbstractArchitecture) where {T} """ function toqubo! end @doc raw""" - toqubo_sense!(model::VirtualQUBOModel, ::AbstractArchitecture) where {T} + toqubo_sense!(model::VirtualModel, ::AbstractArchitecture) where {T} Copies `MOI.ObjectiveSense` from `model.source_model` to `model.target_model`. """ function toqubo_sense! end @doc raw""" - toqubo_variables!(model::VirtualQUBOModel{T}) where {T} + toqubo_variables!(model::VirtualModel{T}) where {T} """ function toqubo_variables! end @doc raw""" - toqubo_variables(model::VirtualQUBOModel{T}) where {T} + toqubo_variables(model::VirtualModel{T}) where {T} """ function toqubo_variable end @doc raw""" - toqubo_objective!(model::VirtualQUBOModel, ::AbstractArchitecture) + toqubo_objective!(model::VirtualModel, ::AbstractArchitecture) """ function toqubo_objective! end @doc raw""" - toqubo_objective(model::VirtualQUBOModel, F::VI, ::AbstractArchitecture) - toqubo_objective(model::VirtualQUBOModel{T}, F::SAF{T}, ::AbstractArchitecture) where {T} - toqubo_objective(model::VirtualQUBOModel{T}, F::SQF{T}, ::AbstractArchitecture) where {T} + toqubo_objective(model::VirtualModel, F::VI, ::AbstractArchitecture) + toqubo_objective(model::VirtualModel{T}, F::SAF{T}, ::AbstractArchitecture) where {T} + toqubo_objective(model::VirtualModel{T}, F::SQF{T}, ::AbstractArchitecture) where {T} """ function toqubo_objective end @doc raw""" - toqubo_constraints!(model::VirtualQUBOModel, ::AbstractArchitecture) + toqubo_constraints!(model::VirtualModel, ::AbstractArchitecture) """ function toqubo_constraints! end @doc raw""" - toqubo_constraint(model::VirtualQUBOModel{T}, f::SAF{T}, s::EQ{T}, ::AbstractArchitecture) where {T} - toqubo_constraint(model::VirtualQUBOModel{T}, f::SAF{T}, s::LT{T}, ::AbstractArchitecture) where {T} - toqubo_constraint(model::VirtualQUBOModel{T}, f::SQF{T}, s::EQ{T}, ::AbstractArchitecture) where {T} - toqubo_constraint(model::VirtualQUBOModel{T}, f::SQF{T}, s::LT{T}, ::AbstractArchitecture) where {T} - toqubo_constraint( - ::VirtualQUBOModel{T}, - ::VI, - ::Union{ - MOI.ZeroOne, - MOI.Integer, - MOI.Interval{T}, - MOI.LessThan{T}, - MOI.GreaterThan{T} - }, - ::AbstractArchitecture - ) where {T} + toqubo_constraint + +Returns the pseudo-boolean function associated to a given constraint from the source model. """ function toqubo_constraint end @doc raw""" - toqubo_parse(::VirtualQUBOModel{T}, ::MOI.AbstractFunction, ::AbstractArchitecture) where {T} + toqubo_parse!( + model::VirtualModel{T}, + g::PBO.PBF{VI,T}, + f::MOI.AbstractFunction, + arch::AbstractArchitectur + ) where {T} -Parses the given function into a PBF. -""" function toqubo_parse end +Parses the given MOI function `f` into PBF `g`. +""" function toqubo_parse! end @doc raw""" - toqubo_penalties!(model::VirtualQUBOModel, ::AbstractArchitecture) + toqubo_penalties!(model::VirtualModel, ::AbstractArchitecture) """ function toqubo_penalties! end @doc raw""" """ function toqubo_penalty end @doc raw""" - toqubo_build!(model::VirtualQUBOModel, ::AbstractArchitecture) -""" function toqubo_build! end \ No newline at end of file + toqubo_build!(model::VirtualModel, ::AbstractArchitecture) +""" function toqubo_build! end + +@doc raw""" + toqubo_empty!(model::VirtualModel, ::AbstractArchitecture) +""" function toqubo_empty! end \ No newline at end of file diff --git a/src/compiler/objective.jl b/src/compiler/objective.jl index 6f9021da..3d34aae2 100644 --- a/src/compiler/objective.jl +++ b/src/compiler/objective.jl @@ -1,4 +1,4 @@ -function toqubo_sense!(model::VirtualQUBOModel, ::AbstractArchitecture) +function toqubo_sense!(model::VirtualModel, ::AbstractArchitecture) if MOI.get(model, MOI.ObjectiveSense()) === MOI.MAX_SENSE MOI.set(model.target_model, MOI.ObjectiveSense(), MOI.MAX_SENSE) else @@ -9,41 +9,11 @@ function toqubo_sense!(model::VirtualQUBOModel, ::AbstractArchitecture) return nothing end -function toqubo_objective!(model::VirtualQUBOModel, arch::AbstractArchitecture) +function toqubo_objective!(model::VirtualModel, arch::AbstractArchitecture) F = MOI.get(model, MOI.ObjectiveFunctionType()) f = MOI.get(model, MOI.ObjectiveFunction{F}()) - copy!(model.f, toqubo_objective(model, f, arch)) + toqubo_parse!(model, model.f, f, arch) return nothing end - -function toqubo_objective( - model::VirtualQUBOModel{T}, - vi::VI, - ::AbstractArchitecture, -) where {T} - f = PBO.PBF{VI,T}() - - for (ω, c) in expansion(MOI.get(model, Source(), vi)) - f[ω] += c - end - - return f -end - -function toqubo_objective( - model::VirtualQUBOModel{T}, - f::SAF{T}, - arch::AbstractArchitecture, -) where {T} - return toqubo_parse(model, f, arch) -end - -function toqubo_objective( - model::VirtualQUBOModel{T}, - f::SQF{T}, - arch::AbstractArchitecture, -) where {T} - return toqubo_parse(model, f, arch) -end \ No newline at end of file diff --git a/src/compiler/parse.jl b/src/compiler/parse.jl index d01c6030..2952b2cd 100644 --- a/src/compiler/parse.jl +++ b/src/compiler/parse.jl @@ -1,40 +1,91 @@ -function toqubo_parse(model::VirtualQUBOModel{T}, f::SAF{T}, ::AbstractArchitecture) where {T} +function toqubo_parse( + model::VirtualModel{T}, + f::MOI.AbstractFunction, + s::MOI.AbstractSet, + arch::AbstractArchitecture, +) where {T} g = PBO.PBF{VI,T}() + toqubo_parse!(model, g, f, s, arch) + + return g +end + +function toqubo_parse!( + model::VirtualModel{T}, + g::PBO.PBF{VI,T}, + vi::VI, + ::AbstractArchitecture, +) where {T} + empty!(g) + + for (ω, c) in expansion(model.source[vi]) + g[ω] += c + end + + return nothing +end + +function toqubo_parse!( + model::VirtualModel{T}, + g::PBO.PBF{VI,T}, + f::SAF{T}, + ::AbstractArchitecture, +) where {T} + empty!(g) + sizehint!(g, length(f.terms) + 1) - + for a in f.terms c = a.coefficient x = a.variable + v = model.source[x] - for (ω, d) in expansion(MOI.get(model, Source(), x)) + for (ω, d) in expansion(v) g[ω] += c * d end end g[nothing] += f.constant - return g + return nothing end -function toqubo_parse(model::VirtualQUBOModel{T}, f::SAF{T}, s::EQ{T}, arch::AbstractArchitecture) where {T} - g = toqubo_parse(model, f, arch) +function toqubo_parse!( + model::VirtualModel{T}, + g::PBO.PBF{VI,T}, + f::SAF{T}, + s::EQ{T}, + arch::AbstractArchitecture, +) where {T} + toqubo_parse!(model, g, f, arch) g[nothing] -= s.value - return g + return nothing end -function toqubo_parse(model::VirtualQUBOModel{T}, f::SAF{T}, s::LT{T}, arch::AbstractArchitecture) where {T} - g = toqubo_parse(model, f, arch) - +function toqubo_parse!( + model::VirtualModel{T}, + g::PBO.PBF{VI,T}, + f::SAF{T}, + s::LT{T}, + arch::AbstractArchitecture, +) where {T} + toqubo_parse!(model, g, f, arch) + g[nothing] -= s.upper - return g + return nothing end -function toqubo_parse(model::VirtualQUBOModel{T}, f::SQF{T}, ::AbstractArchitecture) where {T} - g = PBO.PBF{VI,T}() +function toqubo_parse!( + model::VirtualModel{T}, + g::PBO.PBF{VI,T}, + f::SQF{T}, + ::AbstractArchitecture, +) where {T} + empty!(g) sizehint!(g, length(f.quadratic_terms) + length(f.affine_terms) + 1) @@ -42,6 +93,8 @@ function toqubo_parse(model::VirtualQUBOModel{T}, f::SQF{T}, ::AbstractArchitect c = q.coefficient xi = q.variable_1 xj = q.variable_2 + vi = model.source[xi] + vj = model.source[xj] # MOI convetion is to write ScalarQuadraticFunction as # ½ x' Q x + a x + b @@ -50,8 +103,8 @@ function toqubo_parse(model::VirtualQUBOModel{T}, f::SQF{T}, ::AbstractArchitect c /= 2 end - for (ωi, di) in expansion(MOI.get(model, Source(), xi)) - for (ωj, dj) in expansion(MOI.get(model, Source(), xj)) + for (ωi, di) in expansion(vi) + for (ωj, dj) in expansion(vj) g[union(ωi, ωj)] += c * di * dj end end @@ -60,29 +113,42 @@ function toqubo_parse(model::VirtualQUBOModel{T}, f::SQF{T}, ::AbstractArchitect for a in f.affine_terms c = a.coefficient x = a.variable + v = model.source[x] - for (ω, d) in expansion(MOI.get(model, Source(), x)) + for (ω, d) in expansion(v) g[ω] += c * d end end g[nothing] += f.constant - return g + return nothing end -function toqubo_parse(model::VirtualQUBOModel{T}, f::SQF{T}, s::EQ{T}, arch::AbstractArchitecture) where {T} - g = toqubo_parse(model, f, arch) +function toqubo_parse!( + model::VirtualModel{T}, + g::PBO.PBF{VI,T}, + f::SQF{T}, + s::EQ{T}, + arch::AbstractArchitecture, +) where {T} + toqubo_parse!(model, g, f, arch) g[nothing] -= s.value - return g + return nothing end -function toqubo_parse(model::VirtualQUBOModel{T}, f::SQF{T}, s::LT{T}, arch::AbstractArchitecture) where {T} - g = toqubo_parse(model, f, arch) +function toqubo_parse!( + model::VirtualModel{T}, + g::PBO.PBF{VI,T}, + f::SQF{T}, + s::LT{T}, + arch::AbstractArchitecture, +) where {T} + toqubo_parse!(model, g, f, arch) g[nothing] -= s.upper - return g + return nothing end \ No newline at end of file diff --git a/src/compiler/penalties.jl b/src/compiler/penalties.jl index f844bc35..d7b3225d 100644 --- a/src/compiler/penalties.jl +++ b/src/compiler/penalties.jl @@ -1,5 +1,5 @@ -function toqubo_penalties!(model::VirtualQUBOModel{T}, ::AbstractArchitecture) where {T} - # -*- :: Invert Sign :: -*- # +function toqubo_penalties!(model::VirtualModel{T}, ::AbstractArchitecture) where {T} + # Invert Sign s = MOI.get(model, MOI.ObjectiveSense()) === MOI.MAX_SENSE ? -1 : 1 β = one(T) # TODO: This should be made a parameter too? Yes! diff --git a/src/compiler/validation.jl b/src/compiler/validation.jl deleted file mode 100644 index 4c9505f9..00000000 --- a/src/compiler/validation.jl +++ /dev/null @@ -1,40 +0,0 @@ -function isqubo(model::MOI.ModelLike) - F = MOI.get(model, MOI.ObjectiveFunctionType()) - - if !(F <: Union{SQF,SAF,VI}) - return false - end - - S = MOI.get(model, MOI.ObjectiveSense()) - - if !(S === MOI.MAX_SENSE || S === MOI.MIN_SENSE) - return false - end - - # List of Variables - v = Set{VI}(MOI.get(model, MOI.ListOfVariableIndices())) - - for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent()) - if (F === VI && S === MOI.ZeroOne) - for cᵢ in MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) - vᵢ = MOI.get(model, MOI.ConstraintFunction(), cᵢ) - - # Account for variable as binary - delete!(v, vᵢ) - end - else - # Non VariableIndex-in-ZeroOne Constraint - return false - end - end - - if !isempty(v) - # Some variable is not covered by binary constraints - return false - end - - true -end - -isqubo(::QUBOModel) = true -isqubo(::VirtualQUBOModel) = true \ No newline at end of file diff --git a/src/compiler/variables.jl b/src/compiler/variables.jl index add7acf2..9c1296fa 100644 --- a/src/compiler/variables.jl +++ b/src/compiler/variables.jl @@ -1,4 +1,4 @@ -function toqubo_variables!(model::VirtualQUBOModel{T}, ::AbstractArchitecture) where {T} +function toqubo_variables!(model::VirtualModel{T}, ::AbstractArchitecture) where {T} # Set of all source variables Ω = Vector{VI}(MOI.get(model, MOI.ListOfVariableIndices())) @@ -8,7 +8,7 @@ function toqubo_variables!(model::VirtualQUBOModel{T}, ::AbstractArchitecture) w ℝ = Dict{VI,Tuple{Union{T,Nothing},Union{T,Nothing}}}() for ci in MOI.get(model, MOI.ListOfConstraintIndices{VI,MOI.ZeroOne}()) - # -*- Binary Variable 😄 -*- + # Binary Variable x = MOI.get(model, MOI.ConstraintFunction(), ci) # Add to set @@ -16,7 +16,7 @@ function toqubo_variables!(model::VirtualQUBOModel{T}, ::AbstractArchitecture) w end for ci in MOI.get(model, MOI.ListOfConstraintIndices{VI,MOI.Integer}()) - # -*- Integer Variable 🤔 -*- + # Integer Variable x = MOI.get(model, MOI.ConstraintFunction(), ci) # Add to dict as unbounded @@ -24,17 +24,17 @@ function toqubo_variables!(model::VirtualQUBOModel{T}, ::AbstractArchitecture) w end for x in setdiff(Ω, 𝔹, keys(ℤ)) - # -*- Real Variable 😢 -*- + # Real Variable ℝ[x] = (nothing, nothing) end for ci in MOI.get(model, MOI.ListOfConstraintIndices{VI,MOI.Interval{T}}()) - # -*- Interval 😄 -*- + # Interval x = MOI.get(model, MOI.ConstraintFunction(), ci) - I = MOI.get(model, MOI.ConstraintSet(), ci) + s = MOI.get(model, MOI.ConstraintSet(), ci) - a = I.lower - b = I.upper + a = s.lower + b = s.upper if haskey(ℤ, x) ℤ[x] = (a, b) @@ -44,25 +44,25 @@ function toqubo_variables!(model::VirtualQUBOModel{T}, ::AbstractArchitecture) w end for ci in MOI.get(model, MOI.ListOfConstraintIndices{VI,LT{T}}()) - # -*- Upper Bound 🤔 -*- + # Upper Bound x = MOI.get(model, MOI.ConstraintFunction(), ci) - I = MOI.get(model, MOI.ConstraintSet(), ci) + s = MOI.get(model, MOI.ConstraintSet(), ci) - b = I.upper + b = s.upper if haskey(ℤ, x) ℤ[x] = (first(ℤ[x]), b) - elseif haskey(ℝ, xᵢ) + elseif haskey(ℝ, x) ℝ[x] = (first(ℝ[x]), b) end end for ci in MOI.get(model, MOI.ListOfConstraintIndices{VI,GT{T}}()) - # -*- Lower Bound 🤔 -*- + # Lower Bound x = MOI.get(model, MOI.ConstraintFunction(), ci) - I = MOI.get(model, MOI.ConstraintSet(), ci) + s = MOI.get(model, MOI.ConstraintSet(), ci) - a = I.lower + a = s.lower if haskey(ℤ, x) ℤ[x] = (a, last(ℤ[x])) @@ -71,7 +71,7 @@ function toqubo_variables!(model::VirtualQUBOModel{T}, ::AbstractArchitecture) w end end - # -*- Discretize Real Ones 🤔 -*- + # Discretize Real Ones for (x, (a, b)) in ℝ if isnothing(a) || isnothing(b) error("Unbounded variable $(x) ∈ ℝ") @@ -87,33 +87,40 @@ function toqubo_variables!(model::VirtualQUBOModel{T}, ::AbstractArchitecture) w # N ≥ log₂(1 + |b - a| / 4τ) # # where τ is the (absolute) tolerance - # TODO: Add τ as parameter + # TODO: Add τ as parameter (DONE) + # TODO: Move this comment to the documentation let - τ = MOI.get(model, Tol(), x) - e = MOI.get(model, VariableEncoding(), x) - - encode!(e, model, x, a, b, τ) + e = MOI.get(model, Attributes.VariableEncodingMethod(), x) + n = MOI.get(model, Attributes.VariableEncodingBits(), x) + + if !isnothing(n) + encode!(model, e, x, a, b, n) + else + τ = MOI.get(model, Attributes.VariableEncodingATol(), x) + encode!(model, e, x, a, b, τ) + end end end end - # -*- Discretize Integer Variables 🤔 -*- + # Discretize Integer Variables for (x, (a, b)) in ℤ if isnothing(a) || isnothing(b) error("Unbounded variable $(x) ∈ ℤ") else - let e = MOI.get(model, VariableEncoding(), x) - encode!(e, model, x, a, b) + let + e = MOI.get(model, Attributes.VariableEncodingMethod(), x) + encode!(model, e, x, a, b) end end end - # -*- Mirror Boolean Variables 😄 -*- + # Mirror Boolean Variables for x in 𝔹 - encode!(Mirror(), model, x) + encode!(model, Mirror(), x) end return nothing end -function toqubo_variable(model::VirtualQUBOModel, ::AbstractArchitecture) end \ No newline at end of file +function toqubo_variable(model::VirtualModel, ::AbstractArchitecture) end \ No newline at end of file diff --git a/src/lib/pbo/PBF.jl b/src/lib/pbo/PBF.jl index 8d8d50bb..71cc6287 100644 --- a/src/lib/pbo/PBF.jl +++ b/src/lib/pbo/PBF.jl @@ -1,9 +1,6 @@ import Base: isiterable -using MutableArithmetics -const MA = MutableArithmetics - -@doc raw""" +@doc raw""" """ function _parseterm end @@ -23,55 +20,62 @@ _parseterm(::Type{S}, ::Type{T}, x::Tuple{<:Union{Vector{S},Set{S}},T}) where {S (Set{S}(first(x)), last(x)) @doc raw""" -A Pseudo-Boolean Function ``f \in \mathscr{F}`` over some field ``\mathbb{T}`` takes the form + PseudoBooleanFunction{V,T}(Ω::Dict{Union{Set{V},Nothing},T}) where {V,T} + +A Pseudo-Boolean Function[^Boros2002] ``f \in \mathscr{F}`` over some field ``\mathbb{T}`` takes the form ```math -f(\mathbf{x}) = \sum_{\omega \in \Omega\left[f\right]} c_\omega \prod_{j \in \omega} \mathbb{x}_j +f(\mathbf{x}) = \sum_{\omega \in \Omega\left[f\right]} c_\omega \prod_{j \in \omega} x_j ``` -where each ``\Omega\left[{f}\right]`` is the multi-linear representation of ``f`` as a set of terms. Each term is given by a unique set of indices ``\omega \subseteq \mathbb{S}`` related to some coefficient ``c_\omega \in \mathbb{T}``. We say that ``\omega \in \Omega\left[{f}\right] \iff c_\omega \neq 0``. -Variables ``\mathbf{x}_i`` are indeed boolean, thus ``f : \mathbb{B}^{n} \to \mathbb{T}``. +where each ``\Omega\left[{f}\right]`` is the multi-linear representation of ``f`` as a set of terms. +Each term is given by a unique set of indices ``\omega \subseteq \mathbb{S}`` related to some coefficient ``c_\omega \in \mathbb{T}``. +We say that ``\omega \in \Omega\left[{f}\right] \iff c_\omega \neq 0``. +Variables ``x_j`` are boolean, thus ``f : \mathbb{B}^{n} \to \mathbb{T}``. -## References - * [1] Endre Boros, Peter L. Hammer, Pseudo-Boolean optimization, Discrete Applied Mathematics, 2002 [{doi}](https://doi.org/10.1016/S0166-218X(01)00341-9) +[^Boros2002]: + Endre Boros, Peter L. Hammer, **Pseudo-Boolean optimization**, *Discrete Applied Mathematics*, 2002 [{doi}](https://doi.org/10.1016/S0166-218X(01)00341-9) """ -struct PseudoBooleanFunction{S,T} - Ω::Dict{Set{S},T} +struct PseudoBooleanFunction{V,T} + Ω::Dict{Set{V},T} - function PseudoBooleanFunction{S,T}(Ω::Dict{<:Union{Set{S},Nothing},T}) where {S,T} - new{S,T}( - Dict{Set{S},T}(isnothing(ω) ? Set{S}() : ω => c for (ω, c) in Ω if !iszero(c)), + function PseudoBooleanFunction{V,T}(Ω::Dict{<:Union{Set{V},Nothing},T}) where {V,T} + return new{V,T}( + Dict{Set{V},T}(isnothing(ω) ? Set{V}() : ω => c for (ω, c) in Ω if !iszero(c)), ) end - function PseudoBooleanFunction{S,T}(v::Vector) where {S,T} - Ω = Dict{Set{S},T}() + function PseudoBooleanFunction{V,T}(v::Vector) where {V,T} + Ω = Dict{Set{V},T}() for x in v - ω, a = _parseterm(S, T, x) + ω, a = _parseterm(V, T, x) Ω[ω] = get(Ω, ω, zero(T)) + a end - PseudoBooleanFunction{S,T}(Ω) + return PseudoBooleanFunction{V,T}(Ω) end - function PseudoBooleanFunction{S,T}(x::Base.Generator) where {S,T} - PseudoBooleanFunction{S,T}(collect(x)) + function PseudoBooleanFunction{V,T}(x::Base.Generator) where {V,T} + return PseudoBooleanFunction{V,T}(collect(x)) end - function PseudoBooleanFunction{S,T}(x::Vararg{Any}) where {S,T} - PseudoBooleanFunction{S,T}(collect(x)) + function PseudoBooleanFunction{V,T}(x::Vararg{Any}) where {V,T} + return PseudoBooleanFunction{V,T}(collect(x)) end - function PseudoBooleanFunction{S,T}() where {S,T} - new{S,T}(Dict{Set{S},T}()) + function PseudoBooleanFunction{V,T}() where {V,T} + return new{V,T}(Dict{Set{V},T}()) end end -# -*- Alias -*- -const PBF{S,T} = PseudoBooleanFunction{S,T} +# Alias +const PBF{V,T} = PseudoBooleanFunction{V,T} + +# Broadcast as scalar +Base.broadcastable(f::PBF) = Ref(f) -#-*- Copy -*- +# Copy function Base.copy!(f::PBF{S,T}, g::PBF{S,T}) where {S,T} sizehint!(f, length(g)) copy!(f.Ω, g.Ω) @@ -83,7 +87,7 @@ function Base.copy(f::PBF{S,T}) where {S,T} return copy!(PBF{S,T}(), f) end -# -*- Iterator & Length -*- +# Iterator & Length Base.keys(f::PBF) = keys(f.Ω) Base.values(f::PBF) = values(f.Ω) Base.length(f::PBF) = length(f.Ω) @@ -92,32 +96,38 @@ Base.isempty(f::PBF) = isempty(f.Ω) Base.iterate(f::PBF) = iterate(f.Ω) Base.iterate(f::PBF, i::Integer) = iterate(f.Ω, i) -Base.haskey(f::PBF{S}, k::Set{S}) where {S} = haskey(f.Ω, k) -Base.haskey(f::PBF{S}, k::S) where {S} = haskey(f.Ω, Set{S}([k])) -Base.haskey(f::PBF{S}, ::Nothing) where {S} = haskey(f.Ω, Set{S}()) +Base.haskey(f::PBF{S}, ω::Set{S}) where {S} = haskey(f.Ω, ω) +Base.haskey(f::PBF{S}, ξ::S) where {S} = haskey(f, Set{S}([ξ])) +Base.haskey(f::PBF{S}, ::Nothing) where {S} = haskey(f, Set{S}()) -# -*- Indexing: Get -*- +# Indexing: Get # Base.getindex(f::PBF{S,T}, ω::Set{S}) where {S,T} = get(f.Ω, ω, zero(T)) Base.getindex(f::PBF{S}, η::Vector{S}) where {S} = getindex(f, Set{S}(η)) -Base.getindex(f::PBF{S}, ξ::S...) where {S} = getindex(f, Set{S}(ξ)) +Base.getindex(f::PBF{S}, ξ::S) where {S} = getindex(f, Set{S}([ξ])) Base.getindex(f::PBF{S}, ::Nothing) where {S} = getindex(f, Set{S}()) -# -*- Indexing: Set -*- +# Indexing: Set # function Base.setindex!(f::PBF{S,T}, c::T, ω::Set{S}) where {S,T} if !iszero(c) setindex!(f.Ω, c, ω) - elseif haskey(f.Ω, ω) - delete!(f.Ω, ω) + elseif haskey(f, ω) + delete!(f, ω) end return c end Base.setindex!(f::PBF{S,T}, c::T, η::Vector{S}) where {S,T} = setindex!(f, c, Set{S}(η)) -Base.setindex!(f::PBF{S,T}, c::T, ξ::S...) where {S,T} = setindex!(f, c, Set{S}(ξ)) +Base.setindex!(f::PBF{S,T}, c::T, ξ::S) where {S,T} = setindex!(f, c, Set{S}([ξ])) Base.setindex!(f::PBF{S,T}, c::T, ::Nothing) where {S,T} = setindex!(f, c, Set{S}()) -# -*- Properties -*- +# Indexing: Delete # +Base.delete!(f::PBF{S}, ω::Set{S}) where {S} = delete!(f.Ω, ω) +Base.delete!(f::PBF{S}, η::Vector{S}) where {S} = delete!(f, Set{S}(η)) +Base.delete!(f::PBF{S}, k::S) where {S} = delete!(f, Set{S}([k])) +Base.delete!(f::PBF{S}, ::Nothing) where {S} = delete!(f, Set{S}()) + +# Properties Base.size(f::PBF{S,T}) where {S,T} = (length(f),) function Base.sizehint!(f::PBF, n::Integer) @@ -126,7 +136,7 @@ function Base.sizehint!(f::PBF, n::Integer) return f end -# -*- Comparison: (==, !=, ===, !==) -*- # +# Comparison: (==, !=, ===, !==) # Base.:(==)(f::PBF{S,T}, g::PBF{S,T}) where {S,T} = f.Ω == g.Ω Base.:(==)(f::PBF{S,T}, a::T) where {S,T} = isscalar(f) && (f[nothing] == a) Base.:(!=)(f::PBF{S,T}, g::PBF{S,T}) where {S,T} = f.Ω != g.Ω @@ -151,7 +161,7 @@ Base.one(::Type{PBF{S,T}}) where {S,T} = PBF{S,T}(one(T)) Base.isone(f::PBF) = isscalar(f) && isone(f[nothing]) Base.round(f::PBF{S,T}; kw...) where {S,T} = PBF{S,T}(ω => round(c; kw...) for (ω, c) in f) -# -*- Arithmetic: (+) -*- +# Arithmetic: (+) function Base.:(+)(f::PBF{S,T}, g::PBF{S,T}) where {S,T} h = copy(f) @@ -177,7 +187,7 @@ end Base.:(+)(f::PBF{S,T}, c) where {S,T} = +(f, convert(T, c)) Base.:(+)(c, f::PBF) = +(f, c) -# -*- Arithmetic: (-) -*- +# Arithmetic: (-) function Base.:(-)(f::PBF{S,T}) where {S,T} return PBF{S,T}(Dict{Set{S},T}(ω => -c for (ω, c) in f)) end @@ -217,7 +227,7 @@ end Base.:(-)(c, f::PBF{S,T}) where {S,T} = -(convert(T, c), f) Base.:(-)(f::PBF{S,T}, c) where {S,T} = -(f, convert(T, c)) -# -*- Arithmetic: (*) -*- +# Arithmetic: (*) function Base.:(*)(f::PBF{S,T}, g::PBF{S,T}) where {S,T} h = zero(PBF{S,T}) m = length(f) @@ -265,7 +275,7 @@ end Base.:(*)(f::PBF{S,T}, a) where {S,T} = *(f, convert(T, a)) Base.:(*)(a, f::PBF) = *(f, a) -# -*- Arithmetic: (/) -*- +# Arithmetic: (/) function Base.:(/)(f::PBF{S,T}, a::T) where {S,T} if iszero(a) throw(DivideError()) @@ -276,7 +286,7 @@ end Base.:(/)(f::PBF{S,T}, a) where {S,T} = /(f, convert(T, a)) -# -*- Arithmetic: (^) -*- +# Arithmetic: (^) function Base.:(^)(f::PBF{S,T}, n::Integer) where {S,T} if n < 0 throw(DivideError()) @@ -297,7 +307,7 @@ function Base.:(^)(f::PBF{S,T}, n::Integer) where {S,T} end end -# -*- Arithmetic: Evaluation -*- +# Arithmetic: Evaluation function (f::PBF{S,T})(x::Dict{S,U}) where {S,T,U<:Integer} g = PBF{S,T}() @@ -333,7 +343,7 @@ function (f::PBF{S})() where {S} return f(Dict{S,Int}()) end -# -*- Type conversion -*- +# Type conversion function Base.convert(U::Type{<:T}, f::PBF{<:Any,T}) where {T} if isempty(f) return zero(U) @@ -343,264 +353,3 @@ function Base.convert(U::Type{<:T}, f::PBF{<:Any,T}) where {T} error("Can't convert non-constant Pseudo-boolean Function to scalar type '$U'") end end - - -# -*- MA Mutable Arithmetics -*- - -MA.mutability(::Type{<:PBF{S,T}}) where {S,T} = MA.IsMutable() - - -# -*- MA Mutable Copy: -*- -function MA.mutable_copy(f::PBF{S,T}) where{S,T} - f_copy = PBO.PBF{S,T}() - for (ω,c) in f - f_copy[ω] = MA.copy_if_mutable(f[ω]) - end - return f_copy -end - -# -*- MA isequal_canonical: -*- -function MA.isequal_canonical( - f::PBF{S,T}, - g::PBF{S,T} - ) where {S,T} - - return f == g - -end - - -# -*- MA Arithmetic: zero,one -*- -function MA.operate!( - ::typeof(zero), - f::PBF{S,T} - ) where {S,T} - - empty!(f) - - return f -end - -function MA.operate!( - ::typeof(one), - f::PBF{S,T} - ) where {S,T} - - empty!(f) - f[nothing] = one(T) - - return f -end - - -# -*- MA Arithmetic: (+) -*- - -function MA.operate!( - ::typeof(+), - f::PBF{S,T}, - g::PBF{S,T} - ) where {S,T} - - for(ω, c) in g - f[ω] += c - end - - return f -end - -function MA.operate!( - ::typeof(+), - f::PBF{<:Any,T}, - c::T - ) where {T} - if iszero(c) - return f - end - - f[nothing] += c - - return f -end - -# -*- MA Arithmetic: (-) -*- - - -function MA.operate!( - ::typeof(-), - f::PBF{S,T} - ) where {S,T} - - for(ω, c) in f - f[ω] = -c - end - - return f -end - -function MA.operate!( - ::typeof(-), - f::PBF{S,T}, - g::PBF{S,T} - ) where {S,T} - - for(ω, c) in g - f[ω] -= c - end - - return f -end - - -function MA.operate!( - ::typeof(-), - f::PBF{<:Any,T}, - c::T - ) where {T} - if iszero(c) - return f - end - - f[nothing] -= c - - return f -end - -# -*- MA Arithmetic: (*) -*- - -function MA.operate!( - ::typeof(*), - f::PBF{S,T}, - g::PBF{S,T} - ) where {S,T} - - h = PBF{S,T}() - - if isempty(f) || isempty(g) - return h - else - for (ωᵢ, cᵢ) in f, (ωⱼ, cⱼ) in g - h[union(ωᵢ, ωⱼ)] += cᵢ * cⱼ - end - - return h - end - -end - -function MA.operate!( - ::typeof(*), - f::PBF{S,T}, - a::T - ) where {S,T} - - if iszero(a) - empty!(f) - else - for(ω, c) in f - f[ω] = c*a - end - end - - return f -end - -# -*- MA Arithmetic: (add_mul) -*- -function MA.operate!( - ::typeof(add_mul), - f::PBF{S,T}, - m::T, - g::PBF{S,T} - ) where {S,T} - - for(ω, c) in g - f[ω] += m*c - end - return f - -end - -# -*- MA Arithmetic: (sub_mul) -*- -function MA.operate!( - ::typeof(sub_mul), - f::PBF{S,T}, - m::T, - g::PBF{S,T} - ) where {S,T} - - for(ω, c) in g - f[ω] -= m*c - end - return f - -end - -# -*- MA Arithmetic: (/) -*- - -function MA.operate!( - ::typeof(/), - f::PBF{S,T}, - a::T - ) where {S,T} - - if iszero(a) - throw(DivideError()) - else - for(ω, c) in f - f[ω] = c / a - end - end - - return f -end - -# -*- MA Arithmetic: (^) -*- -function MA.operate!( - ::typeof(^), - f::PBF{S,T}, - n::Integer - ) where {S,T} - - if n < 0 - throw(DivideError()) - elseif n == 0 - return MA.operate!(one, f) - elseif n == 1 - return f - elseif n == 2 - for (ω,c) in f - f[ω] = c*c - end - else - for (ω,c) in f - f[ω] = c*c - end - - if iseven(n) - return f^(n ÷ 2) - else - return f^(n ÷ 2) * f - end - end -end - - -function MA.operate_to!( - output::PBF{S,T}, - op::Union{typeof(+),typeof(-)}, - f::PBF{S,T}, - g::PBF{S,T}, -) where {S,T} - empty!(output) - MA.operate!(+, output, f) - MA.operate!(op, output, g) - return output -end - - - -# -*- MA Arithmetic: Evaluation -*- - - - -# -*- MA Type conversion -*- - diff --git a/src/lib/pbo/PBO.jl b/src/lib/pbo/PBO.jl index 895bc6b7..e39470c0 100644 --- a/src/lib/pbo/PBO.jl +++ b/src/lib/pbo/PBO.jl @@ -2,12 +2,8 @@ module PBO using Random -# -*- Imports -*- -import QUBOTools: varcmp, qubo, variables, variable_map, variable_set, variable_inv - -# -*- Exports -*- -# export PBF -# export ×, ∂, Δ, Θ, δ, ϵ, Ω +# Imports +import QUBOTools: varlt, qubo, variables, variable_map, variable_set, variable_inv include("PBF.jl") include("tools.jl") diff --git a/src/lib/pbo/interface.jl b/src/lib/pbo/interface.jl index 7b5d47d2..c1586c81 100644 --- a/src/lib/pbo/interface.jl +++ b/src/lib/pbo/interface.jl @@ -1,9 +1,56 @@ + +@doc raw""" + gap(f::PBF{S, T}; bound::Symbol=:loose) where {S, T} + +Computes the least upper bound for the greatest variantion possible under some `` f \in \mathscr{F} `` i. e. + +```math +\begin{array}{r l} + \min & M \\ + \text{s.t.} & \left|{f(\mathbf{x}) - f(\mathbf{y})}\right| \le M ~~ \forall \mathbf{x}, \mathbf{y} \in \mathbb{B}^{n} +\end{array} +``` + +A simple approach, avaiable using the `bound=:loose` parameter, is to define +```math +M \triangleq \sum_{\omega \neq \varnothing} \left|{c_\omega}\right| +``` +""" +function gap end + +const δ = gap + +@doc raw""" + sharpness(f::PBF{S, T}; bound::Symbol=:loose, tol::T = T(1e-6)) where {S, T} +""" +function sharpness end + +const ϵ = sharpness + +@doc raw""" + derivative(f::PBF{V,T}, x::V) where {V,T} + +The partial derivate of function ``f \in \mathscr{F}`` with respect to the ``x`` variable. + +```math + \Delta_i f(\mathbf{x}) = \frac{\partial f(\mathbf{x})}{\partial \mathbf{x}_i} = + \sum_{\omega \in \Omega\left[{f}\right] \setminus \left\{{i}\right\}} + c_{\omega \cup \left\{{i}\right\}} \prod_{k \in \omega} \mathbf{x}_k +``` +""" +function derivative end + +const Δ = derivative +const ∂ = derivative + @doc raw""" gradient(f::PBF) Computes the gradient of ``f \in \mathscr{F}`` where the ``i``-th derivative is given by [`derivative`](@ref). """ function gradient end +const ∇ = gradient + @doc raw""" residual(f::PBF{S, T}, x::S) where {S, T} @@ -19,7 +66,7 @@ The residual of ``f \in \mathscr{F}`` with respect to ``x``. const Θ = residual @doc raw""" - discretize(f::PBF{S, T}; tol::T) where {S, T} + discretize(f::PBF{V,T}; tol::T) where {V,T} For a given function ``f \in \mathscr{F}`` written as @@ -35,4 +82,10 @@ computes an approximate function ``g : \mathbb{B}^{n} \to \mathbb{Z}`` such tha This is done by rationalizing every coefficient ``c_\omega`` according to some tolerance `tol`. -""" function discretize end \ No newline at end of file +""" function discretize end + +@doc raw""" + discretize!(f::PBF{V,T}; tol::T) where {V,T} + +In-place version of [`discretize`](@ref). +""" function discretize! end \ No newline at end of file diff --git a/src/lib/pbo/print.jl b/src/lib/pbo/print.jl index 0357182f..2bcbc365 100644 --- a/src/lib/pbo/print.jl +++ b/src/lib/pbo/print.jl @@ -3,7 +3,7 @@ function showvar end showvar(x::Any) = x -showvar(s::Set) = showvar.(sort(collect(s); lt=PBO.varcmp)) +showvar(s::Set) = showvar.(sort(collect(s); lt=PBO.varlt)) showvar(x::Integer, v::Symbol = :x) = join([v; Char(0x2080) .+ reverse(digits(x))]) function showterm(ω::Set{S}, c::T, isfirst::Bool) where {S, T} @@ -19,11 +19,13 @@ function showterm(ω::Set{S}, c::T, isfirst::Bool) where {S, T} end function Base.show(io::IO, f::PBF{<:Any, T}) where T + Ω = sort!(collect(f); lt=(x,y) -> varlt(first(x), first(y))) + print(io, if isempty(f) zero(T) else - join(showterm(ω, c, isone(i)) for (i, (ω, c)) in enumerate(f)) + join(showterm(ω, c, isone(i)) for (i, (ω, c)) in enumerate(Ω)) end ) end \ No newline at end of file diff --git a/src/lib/pbo/quadratization.jl b/src/lib/pbo/quadratization.jl index 65e85662..ab14e285 100644 --- a/src/lib/pbo/quadratization.jl +++ b/src/lib/pbo/quadratization.jl @@ -1,205 +1,178 @@ -# -*- :: Quadratization :: -*- - +# :: Quadratization :: # abstract type QuadratizationMethod end -_aux(::Type{<:QuadratizationMethod}, ::Int) = 0 -_nst(::Type{<:QuadratizationMethod}, ::Int) = 0 - -struct Quadratization{T<:QuadratizationMethod} - deg::Int # Initial Degree - aux::Int # Auxiliary variables - nst::Int # Non-Submodular Terms +struct Quadratization{Q<:QuadratizationMethod} + stable::Bool - function Quadratization{T}(deg::Int) where {T<:QuadratizationMethod} - return new{T}( - deg, - _aux(T, deg), - _nst(T, deg), - ) + function Quadratization{Q}(stable::Bool = false) where {Q<:QuadratizationMethod} + return new{Q}(stable) end end @doc raw""" - @quadratization(name, aux, nst) - -Defines a new quadratization technique. -""" -macro quadratization(name, aux, nst) - aux_func = if aux isa Integer - quote _aux(::Type{$(esc(name))}, ::Int) = $(aux) end - else - quote _aux(::Type{$(esc(name))}, deg::Int) = $(aux)(deg) end - end - - nst_func = if nst isa Integer - quote _nst(::Type{$(esc(name))}, ::Int) = $(nst) end - else - quote _nst(::Type{$(esc(name))}, deg::Int) = $(nst)(deg) end - end + quadratize!(aux::Function, f::PBF{S, T}, ::Quadratization{Q}) where {S,T,Q} - return quote - struct $(esc(name)) <: QuadratizationMethod end - - # Ancillary Variables - $(aux_func) - - # Non-Submodular Terms - $(nst_func) - end -end - -@doc raw""" - infer_quadratization(f::PBF) -""" -function infer_quadratization(f::PBF) - k = degree(f) - - if k <= 2 - return nothing - else - # Without any extra knowledge, it is better to - # quadratize term-by-term - return Quadratization{TBT}(k) - end -end - -@doc raw""" - quadratize(aux::Function, f::PBF{S, T}, ::Quadratization) where {S, T} - -Quadratizes a given PBF, i.e. creates a function ``g \in \mathscr{F}^{2}`` from ``f \in \mathscr{F}^{k}, k \ge 3``. +Quadratizes a given PBF in-place, i.e. applies a mapping ``Q : \mathscr{F}^{k} \to \mathscr{F}^{2}``. ```julia aux(::Nothing)::S aux(::Integer)::Vector{S} ``` - -## Submodularity - -A function ``f : 2^{S} \to \mathbb{R}`` is said to be submodular if -```math -f(X \cup Y) + f(X \cap Y) \le f(X) + f(Y) \forall X, Y \subset S -``` -""" function quadratize end +""" function quadratize! end @doc raw""" - Quadratization{TBT}(::Int) + Quadratization{NTR_KZFD}(stable::Bool = false) -Term-by-term quadratization. Employs other inner methods. -""" +NTR-KZFD (Kolmogorov & Zabih, 2004; Freedman & Drineas, 2005) -@quadratization(TBT, 0, 0) +!!! info + Introduces one new variable and no non-submodular terms. +""" struct NTR_KZFD <: QuadratizationMethod end + +function quadratize!( + aux::Function, + f::PBF{S,T}, + ω::Set{S}, + c::T, + ::Quadratization{NTR_KZFD}, +) where {S,T} + # Degree + k = length(ω) -@doc raw""" - Quadratization{NTR_KZFD}(::Int) + # Fast-track + k < 3 && return nothing -NTR-KZFD (Kolmogorov & Zabih, 2004; Freedman & Drineas, 2005) -""" + # Variables + s = aux()::S -@quadratization(NTR_KZFD, 1, 0) + # Stabilize + # NOTE: This method is stable by construction -function quadratize(aux::Function, ω::Set{S}, c::T, ::Quadratization{NTR_KZFD}) where {S, T} - # -* Degree *- - k = length(ω) - s = aux()::S + # Quadratization + delete!(f, ω) - g = PBF{S,T}(s => -c * (k - 1)) + f[s] += -c * (k - 1) for i ∈ ω - η = Set{S}([s, i]) - - g[η] += c + f[i×s] += c end - return g + return nothing end @doc raw""" - Quadratization{PTR_BG}(::Int) + Quadratization{PTR_BG}(stable::Bool = false) PTR-BG (Boros & Gruber, 2014) -""" - -@quadratization( - PTR_BG, - k -> k - 2, - k -> k - 1, -) -function quadratize(aux::Function, ω::Set{S}, c::T, ::Quadratization{PTR_BG}) where {S, T} - # -* Degree *- +!!! info + Introduces ``k - 2`` new variables and ``k - 1`` non-submodular terms. +""" struct PTR_BG <: QuadratizationMethod end + +function quadratize!( + aux::Function, + f::PBF{S,T}, + ω::Set{S}, + c::T, + quad::Quadratization{PTR_BG}, +) where {S,T} + # Degree k = length(ω) - # -* Variables *- + # Fast-track + k < 3 && return nothing + + # Variables s = aux(k - 2)::Vector{S} - b = sort(collect(ω); lt = varcmp)::Vector{S} + b = collect(ω)::Vector{S} + + # Stabilize + quad.stable && sort!(b; lt = varlt) + + # Quadratization + delete!(f, ω) - # -*- Quadratization -*- - f = PBF{S, T}(b[k] × b[k - 1] => c) + f[b[k]×b[k-1]] += c - for i = 1:(k - 2) + for i = 1:(k-2) f[s[i]] += c * (k - i - 1) - f[s[i] × b[i]] += c + f[s[i]×b[i]] += c - for j = (i + 1):k - f[s[i] × b[j]] -= c + for j = (i+1):k + f[s[i]×b[j]] -= c end - end + end - return f + return nothing end -function quadratize(aux::Function, ω::Set{S}, c::T) where {S, T} +@doc raw""" + Quadratization{TERM_BY_TERM}(stable::Bool = false) + +Term-by-term quadratization. Employs other inner methods, specially [`NTR_KZFD`](@ref) and [`PTR_BG`](@ref). +""" struct TERM_BY_TERM <: QuadratizationMethod end + +function quadratize!( + aux::Function, + f::PBF{S,T}, + ω::Set{S}, + c::T, + quad::Quadratization{TERM_BY_TERM}, +) where {S,T} if c < zero(T) - return quadratize( - aux, - ω, - c, - Quadratization{NTR_KZFD}(length(ω)), - ) + quadratize!(aux, f, ω, c, Quadratization{NTR_KZFD}(quad.stable)) else - return quadratize( - aux, - ω, - c, - Quadratization{PTR_BG}(length(ω)), - ) + quadratize!(aux, f, ω, c, Quadratization{PTR_BG}(quad.stable)) end -end - -function quadratize(aux::Function, f::PBF{S, T}, ::Quadratization{TBT}) where {S, T} - g = PBF{S, T}() - sizehint!(g, length(f)) + return nothing +end - for (ω, c) in f - k = length(ω) +function quadratize!( + aux::Function, + f::PBF{S,T}, + quad::Quadratization{TERM_BY_TERM}, +) where {S,T} + # Collect Terms + Ω = collect(f) - if k <= 2 - g[ω] += c - else - h = quadratize( - aux, - ω, - c, - Quadratization{PTR_BG}(k), - ) + # Stable Quadratization + quad.stable && sort!(Ω; by = first, lt = varlt) - for (η, a) in h - g[η] += a - end - end + for (ω, c) in Ω + quadratize!(aux, f, ω, c, Quadratization{PTR_BG}(quad.stable)) end - return g + return nothing end -function quadratize(aux::Function, f::PBF) - quad = infer_quadratization(f) - if isnothing(quad) - return f +@doc raw""" + Quadratization{INFER}(stable::Bool = false) +""" struct INFER <: QuadratizationMethod end + +@doc raw""" + infer_quadratization(f::PBF) +""" +function infer_quadratization(f::PBF, stable::Bool = false) + k = degree(f) + + if k <= 2 + return nothing else - return quadratize(aux, f, quad) + # Without any extra knowledge, it is better to + # quadratize term-by-term + return Quadratization{TERM_BY_TERM}(stable) end -end \ No newline at end of file +end + +function quadratize!(aux::Function, f::PBF, quad::Quadratization{INFER}) + quadratize!(aux, f, infer_quadratization(f, quad.stable)) + + return nothing +end + +function quadratize!(::Function, ::PBF, ::Nothing) + return nothing +end diff --git a/src/lib/pbo/tools.jl b/src/lib/pbo/tools.jl index 28f1ac79..ef1d7470 100644 --- a/src/lib/pbo/tools.jl +++ b/src/lib/pbo/tools.jl @@ -1,6 +1,6 @@ -# -*- Relaxed Greatest Common Divisor -*- +# Relaxed Greatest Common Divisor @doc raw""" - relaxed_gcd(x::T, y::T; tol::T = T(1e-6)) where {T <: AbstractFloat} + relaxed_gcd(x::T, y::T; tol::T = T(1e-6)) where {T} We define two real numbers ``x`` and ``y`` to be ``\tau``-comensurable if, for some ``\tau > 0`` there exists a continued fractions convergent ``p_{k} \div q_{k}`` such that @@ -8,7 +8,7 @@ We define two real numbers ``x`` and ``y`` to be ``\tau``-comensurable if, for s \left| {q_{k} x - p_{k} y} \right| \le \tau ``` """ -function relaxed_gcd(x::T, y::T; tol::T = 1e-6) where {T<:Number} +function relaxed_gcd(x::T, y::T; tol::T = 1e-6) where {T} x_ = abs(x) y_ = abs(y) @@ -23,132 +23,102 @@ function relaxed_gcd(x::T, y::T; tol::T = 1e-6) where {T<:Number} end end -function relaxed_gcd(a::AbstractArray{T}; tol::T = 1e-6) where {T<:Number} +function relaxed_gcd(a::AbstractArray{T}; tol::T = 1e-6) where {T} if length(a) == 0 - one(T) + return one(T) elseif length(a) == 1 - first(a) + return first(a) else - reduce((x, y) -> relaxed_gcd(x, y; tol = tol), a) + return reduce((x, y) -> relaxed_gcd(x, y; tol = tol), a) end end -# -*- Variable Terms -*- -varmul(x::S, y::S) where {S} = Set{S}([x, y]) -varmul(x::Set{S}, y::S) where {S} = push!(copy(x), y) -varmul(x::S, y::Set{S}) where {S} = push!(copy(y), x) -varmul(x::Set{S}, y::Set{S}) where {S} = union(x, y) +# Variable Terms +varmul(x::V, y::V) where {V} = Set{V}([x, y]) +varmul(x::Set{V}, y::V) where {V} = push!(copy(x), y) +varmul(x::V, y::Set{V}) where {V} = push!(copy(y), x) +varmul(x::Set{V}, y::Set{V}) where {V} = union(x, y) -const × = varmul # \times -const ≺ = varcmp # \prec, from QUBOTools +const × = varmul # \times[tab] +const ≺ = varlt # \prec[tab] @doc raw""" """ -function degree end # TODO: memoize +function degree end degree(f::PBF) = maximum(length.(keys(f)); init = 0) -# -*- Gap & Penalties -*- +# Gap & Penalties @doc raw""" """ function lowerbound end -lowerbound(f::PBF; bound::Symbol = :loose) = lowerbound(f, Val(bound)) -lowerbound(f::PBF{<:Any,T}, ::Val{:loose}) where {T} = - sum(c < zero(T) || isempty(ω) ? c : zero(T) for (ω, c) in f) +function lowerbound(f::PBF; bound::Symbol = :loose) + return lowerbound(f, Val(bound)) +end + +function lowerbound(f::PBF{<:Any,T}, ::Val{:loose}) where {T} + return sum(c < zero(T) || isempty(ω) ? c : zero(T) for (ω, c) in f) +end @doc raw""" """ function upperbound end -upperbound(f::PBF; bound::Symbol = :loose) = upperbound(f, Val(bound)) -upperbound(f::PBF{<:Any,T}, ::Val{:loose}) where {T} = - sum(c > zero(T) || isempty(ω) ? c : zero(T) for (ω, c) in f) - -@doc raw""" -""" function bounds end +function upperbound(f::PBF; bound::Symbol = :loose) + return upperbound(f, Val(bound)) +end -bounds(f::PBF; bound::Symbol = :loose) = - (lowerbound(f; bound = bound), upperbound(f; bound = bound)) +function upperbound(f::PBF{<:Any,T}, ::Val{:loose}) where {T} + return sum(c > zero(T) || isempty(ω) ? c : zero(T) for (ω, c) in f) +end @doc raw""" - gap(f::PBF{S, T}; bound::Symbol=:loose) where {S, T} - -Computes the least upper bound for the greatest variantion possible under some `` f \in \mathscr{F} `` i. e. - -```math -\begin{array}{r l} - \min & M \\ - \text{s.t.} & \left|{f(\mathbf{x}) - f(\mathbf{y})}\right| \le M ~~ \forall \mathbf{x}, \mathbf{y} \in \mathbb{B}^{n} -\end{array} -``` +""" function bounds end -A simple approach, avaiable using the `bound=:loose` parameter, is to define -```math -M \triangleq \sum_{\omega \neq \varnothing} \left|{c_\omega}\right| -``` -""" -function gap end +function bounds(f::PBF; bound::Symbol = :loose) + return (lowerbound(f; bound), upperbound(f; bound)) +end -gap(f::PBF; bound::Symbol = :loose) = gap(f, Val(bound)) -gap(f::PBF{<:Any,T}, ::Val{:loose}) where {T} = - sum(abs(c) for (ω, c) in f if !isempty(ω); init = zero(T)) -gap(::PBF, ::Val{:tight}) = error("Not Implemented: See [1] sec 5.1.1 Majorization") +function gap(f::PBF; bound::Symbol = :loose) + return gap(f, Val(bound)) +end -const δ = gap +function gap(f::PBF{V,T}, ::Val{:loose}) where {V,T} + return sum(abs(c) for (ω, c) in f if !isempty(ω); init = zero(T)) +end -@doc raw""" - sharpness(f::PBF{S, T}; bound::Symbol=:loose, tol::T = T(1e-6)) where {S, T} -""" -function sharpness end +function gap(::PBF, ::Val{:tight}) + error("Not Implemented: See [1] sec 5.1.1 Majorization") +end -function sharpness(f::PBF{S,T}; bound::Symbol = :loose, tol::T = 1e-6) where {S,T} +function sharpness(f::PBF{V,T}; bound::Symbol = :loose, tol::T = 1e-6) where {V,T} return sharpness(f, Val(bound), tol) end -function sharpness(::PBF{S,T}, ::Val{:none}, ::T) where {S,T} +function sharpness(::PBF{V,T}, ::Val{:none}, ::T) where {V,T} return one(T) end -function sharpness(f::PBF{S,T}, ::Val{:loose}, tol::T = 1E-6) where {S,T} +function sharpness(f::PBF{V,T}, ::Val{:loose}, tol::T = 1E-6) where {V,T} return relaxed_gcd(collect(values(f)); tol = tol)::T end -const ϵ = sharpness - -@doc raw""" - derivative(f::PBF{S, T}, i::S) where {S, T} - derivative(f::PBF{S, T}, i::Int) where {S, T} - -The partial derivate of function ``f \in \mathscr{F}`` with respect to the ``i``-th variable. - -```math - \Delta_i f(\mathbf{x}) = \frac{\partial f(\mathbf{x})}{\partial \mathbf{x}_i} = - \sum_{\omega \in \Omega\left[{f}\right] \setminus \left\{{i}\right\}} - c_{\omega \cup \left\{{i}\right\}} \prod_{k \in \omega} \mathbf{x}_k -``` -""" -function derivative end - -function derivative(f::PBF{S,T}, s::S) where {S,T} - return PBF{S,T}(ω => f[ω×s] for ω ∈ keys(f) if (s ∉ ω)) +function derivative(f::PBF{V,T}, x::V) where {V,T} + return PBF{V,T}(ω => f[ω×x] for ω ∈ keys(f) if (x ∉ ω)) end -const Δ = derivative -const ∂ = derivative - -function gradient(f::PBF{S,<:Any}, x::Vector{S}) where {S} - return [derivative(f, x) for (s, _) in varmap(f)] +function gradient(f::PBF{V}, x::Vector{V}) where {V} + return derivative.(f, x) end -const ∇ = gradient - -residual(f::PBF{S,T}, i::S) where {S,T} = PBF{S,T}(ω => c for (ω, c) ∈ keys(f) if (i ∉ ω)) -residual(f::PBF, i::Int) = residual(f, varinv(f)[i]) +function residual(f::PBF{V,T}, x::V) where {V,T} + return PBF{V,T}(ω => c for (ω, c) ∈ keys(f) if (x ∉ ω)) +end -function discretize(f::PBF{S,T}; tol::T = 1E-6) where {S,T} +function discretize(f::PBF{V,T}; tol::T = 1E-6) where {V,T} return discretize!(copy(f); tol = tol) end -function discretize!(f::PBF{S,T}; tol::T = 1E-6) where {S,T} +function discretize!(f::PBF{V,T}; tol::T = 1E-6) where {V,T} ε = sharpness(f; bound = :loose, tol = tol) for (ω, c) in f diff --git a/src/lib/pbo/wrapper.jl b/src/lib/pbo/wrapper.jl index cea9e108..b953af85 100644 --- a/src/lib/pbo/wrapper.jl +++ b/src/lib/pbo/wrapper.jl @@ -9,7 +9,7 @@ function qubo(f::PBF{S,T}, ::Type{Dict}) where {S,T} sizehint!(Q, length(f)) for (ω, a) in f - η = sort([x[i] for i ∈ ω]; lt = varcmp) + η = sort([x[i] for i ∈ ω]; lt = varlt) k = length(η) if k == 0 @@ -39,7 +39,7 @@ function qubo(f::PBF{S,T}, ::Type{Matrix}) where {S,T} β = zero(T) for (ω, a) ∈ f - η = sort([x[i] for i ∈ ω]; lt = varcmp) + η = sort([x[i] for i ∈ ω]; lt = varlt) k = length(η) if k == 0 β += a @@ -64,4 +64,4 @@ end variable_map(f::PBF{S}) where {S} = Dict{S,Int}(v => i for (i, v) in enumerate(variables(f))) variable_inv(f::PBF{S}) where {S} = Dict{Int,S}(i => v for (i, v) in enumerate(variables(f))) variable_set(f::PBF{S}) where {S} = reduce(union!, keys(f); init=Set{S}()) -variables(f::PBF) = sort(collect(variable_set(f)); lt = varcmp) +variables(f::PBF) = sort(collect(variable_set(f)); lt = varlt) diff --git a/src/model/attributes.jl b/src/model/attributes.jl deleted file mode 100644 index 91c0310f..00000000 --- a/src/model/attributes.jl +++ /dev/null @@ -1,67 +0,0 @@ -struct Tol <: MOI.AbstractOptimizerAttribute end - -function MOI.get(model::VirtualQUBOModel{T}, ::Tol) where T - return model.settings.atol[nothing]::T -end - -function MOI.set(model::VirtualQUBOModel{T}, ::Tol, atol::T) where T - @assert atol > zero(T) - - model.settings.atol[nothing] = atol -end - -struct Penalty <: MOI.AbstractOptimizerAttribute end - -function MOI.get(model::VirtualQUBOModel{T}, ::Penalty, ci::CI) where T - return model.ρ[ci] -end - -function MOI.set(model::VirtualQUBOModel{T}, ::Penalty, ci::CI, ρ::T) where T - model.ρ[ci] = ρ -end - -function MOI.get(model::VirtualQUBOModel{T}, ::Penalty, vi::VI) where T - return model.ρ[vi] -end - -function MOI.set(model::VirtualQUBOModel{T}, ::Penalty, vi::VI, ρ::T) where T - model.ρ[vi] = ρ -end - -struct GlobalEncoding <: MOI.AbstractOptimizerAttribute end - -function MOI.get(model::VirtualQUBOModel, ::GlobalEncoding) - return model.settings.encoding[nothing] -end - -function MOI.set(model::VirtualQUBOModel, ::GlobalEncoding, e::Encoding) - model.settings.encoding[nothing] = e -end - -struct VariableEncoding <: MOI.AbstractVariableAttribute end - -function MOI.get(model::VirtualQUBOModel, ::VariableEncoding, vi::VI) - if haskey(model.settings.encoding, vi) - return model.settings.encoding[vi] - else - return model.settings.encoding[nothing] - end -end - -function MOI.set(model::VirtualQUBOModel, ::VariableEncoding, vi::VI, e::Encoding) - model.settings.encoding[vi] = e -end - -struct ConstraintEncoding <: MOI.AbstractConstraintAttribute end - -function MOI.get(model::VirtualQUBOModel, ::ConstraintEncoding, ci::CI) - if haskey(model.settings.encoding, ci) - return model.settings.encoding[ci] - else - return MOI.get(model, GlobalEncoding()) - end -end - -function MOI.set(model::VirtualQUBOModel, ::ConstraintEncoding, ci::CI, e::Encoding) - model.settings.encoding[ci] = e -end \ No newline at end of file diff --git a/src/model/prequbo.jl b/src/model/prequbo.jl index a2e174a6..d10ed1e4 100644 --- a/src/model/prequbo.jl +++ b/src/model/prequbo.jl @@ -1,4 +1,3 @@ -# -*- Model: PreQUBOModel -*- # MOIU.@model(PreQUBOModel, # Name of model (MOI.Integer, MOI.ZeroOne), # untyped scalar sets (EQ, LT, GT), # typed scalar sets @@ -11,6 +10,6 @@ MOIU.@model(PreQUBOModel, # Name of model false, # is optimizer? ) -# :: Drop Automatic Constraint Support :: # +# Drop Generic Constraint Support MOI.supports_constraint(::PreQUBOModel{T}, ::Type{SAF{T}}, ::Type{GT{T}}) where {T} = false MOI.supports_constraint(::PreQUBOModel{T}, ::Type{SQF{T}}, ::Type{GT{T}}) where {T} = false \ No newline at end of file diff --git a/src/model/qubo.jl b/src/model/qubo.jl index 101ad886..926a0f4e 100644 --- a/src/model/qubo.jl +++ b/src/model/qubo.jl @@ -1,120 +1,97 @@ -const QUBO_MODEL_BACKEND{T} = QUBOTools.StandardQUBOModel{QUBOTools.BoolDomain,VI,T,Int} - mutable struct QUBOModel{T} <: MOI.ModelLike - model::QUBO_MODEL_BACKEND{T} - sense::MOI.OptimizationSense + objective_function::SQF{T} + objective_sense::MOI.OptimizationSense + variables::Vector{VI} - QUBOModel{T}(args...; kw...) where {T} = new{T}(QUBO_MODEL_BACKEND{T}(), MOI.MIN_SENSE) - QUBOModel(args...; kw...) = QUBOTools{Float64}(args...; kw...) + function QUBOModel{T}() where {T} + return new{T}(SQF{T}(SQT{T}[], SAT{T}[], zero(T)), MOI.MIN_SENSE, VI[]) + end end -QUBOTools.backend(model::QUBOModel) = model.model - -# -*- :: MOI :: -*- # +# MOI Wrapper function MOI.add_variable(model::QUBOModel) - v, _ = MOI.add_constrained_variable(model, MOI.ZeroOne()) + vi, _ = MOI.add_constrained_variable(model, MOI.ZeroOne()) - return v + return vi end function MOI.add_constrained_variable(model::QUBOModel, ::MOI.ZeroOne) - i = QUBOTools.domain_size(model) + 1 - v = MOI.VariableIndex(i) - c = MOI.ConstraintIndex{VI,MOI.ZeroOne}(v.value) + i = MOI.get(model, MOI.NumberOfVariables()) + 1 + vi = MOI.VariableIndex(i) + ci = MOI.ConstraintIndex{VI,MOI.ZeroOne}(i) - push!(QUBOTools.variable_map(model), v => i) - push!(QUBOTools.variable_inv(model), i => v) + push!(model.variables, VI(i)) - return (v, c) + return (vi, ci) end -MOI.is_empty(model::QUBOModel) = isempty(model.model) +function MOI.is_empty(model::QUBOModel) + return isempty(model.variables) && + isempty(model.objective_function.quadratic_terms) && + isempty(model.objective_function.affine_terms) && + iszero(model.objective_function.constant) +end + +function MOI.empty!(model::QUBOModel{T}) where {T} + model.objective_function = SQF{T}(SQT{T}[], SAT{T}[], zero(T)) + model.objective_sense = MOI.MIN_SENSE -function MOI.empty!(model::QUBOModel) - empty!(model.model) + empty!(model.variables) return nothing end -# -*- :: Support :: -*- # +# Support MOI.supports(::QUBOModel{T}, ::MOI.ObjectiveFunction{SQF{T}}) where {T} = true MOI.supports_constraint(::QUBOModel{T}, ::Type{VI}, ::Type{MOI.ZeroOne}) where {T} = true MOI.supports_add_constrained_variable(::QUBOModel{T}, ::Type{MOI.ZeroOne}) where {T} = true function MOI.add_constraint(::QUBOModel, vi::VI, ::MOI.ZeroOne) - MOI.ConstraintIndex{VI,MOI.ZeroOne}(vi.value) + return MOI.ConstraintIndex{VI,MOI.ZeroOne}(vi.value) end -# -*- :: get + set :: -*- # +# get & set function MOI.get(model::QUBOModel, ::MOI.ObjectiveSense) - return model.sense + return model.objective_sense end -function MOI.set(model::QUBOModel, ::MOI.ObjectiveSense, os::MOI.OptimizationSense) - model.sense = os -end - -function MOI.get( - model::QUBOModel{T}, - ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{T}}, -) where {T} - b = QUBOTools.offset(model) - a = MOI.ScalarAffineTerm{T}[] - Q = MOI.ScalarQuadraticTerm{T}[] - - for (i, c) in QUBOTools.linear_terms(model) - xi = QUBOTools.variable_inv(model, i) - - push!(a, MOI.ScalarAffineTerm{T}(c, xi)) - end - - for ((i, j), c) in QUBOTools.quadratic_terms(model) - xi = QUBOTools.variable_inv(model, i) - xj = QUBOTools.variable_inv(model, j) +function MOI.set( + model::QUBOModel, + ::MOI.ObjectiveSense, + objective_sense::MOI.OptimizationSense, +) + model.objective_sense = objective_sense - push!(Q, MOI.ScalarQuadraticTerm{T}(c, xi, xj)) - end + return nothing +end - return MOI.ScalarQuadraticFunction{T}(Q, a, b) +function MOI.get(model::QUBOModel{T}, ::MOI.ObjectiveFunction{SQF{T}}) where {T} + return model.objective_function end -function MOI.set( - model::QUBOModel{T}, - ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{T}}, - f::MOI.ScalarQuadraticFunction{T}, -) where {T} - linear_terms = Dict{VI,T}() - quadratic_terms = Dict{Tuple{VI,VI},T}() - - for a in f.affine_terms - c = a.coefficient - x = a.variable - - linear_terms[x] = get(linear_terms, x, zero(T)) + c - end +function MOI.set(model::QUBOModel{T}, ::MOI.ObjectiveFunction{VI}, vi::VI) where {T} + model.objective_function = SQF{T}(SQT{T}[], SAT{T}[SAT{T}(one(T), vi)], zero(T)) - for q in f.quadratic_terms - c = q.coefficient - xi = q.variable_1 - xj = q.variable_2 + return nothing +end - if xi == xj - linear_terms[xi] = get(linear_terms, xi, zero(T)) + c / 2 - else - quadratic_terms[(xi, xj)] = get(quadratic_terms, (xi, xj), zero(T)) + c - end - end +function MOI.set(model::QUBOModel{T}, ::MOI.ObjectiveFunction{SAF{T}}, f::SAF{T}) where {T} + model.objective_function = SQF{T}(SQT{T}[], copy(f.terms), f.constant) - offset = f.constant + return nothing +end - model.model = QUBO_MODEL_BACKEND{T}(linear_terms, quadratic_terms; offset = offset) +function MOI.set(model::QUBOModel{T}, ::MOI.ObjectiveFunction{SQF{T}}, f::SQF{T}) where {T} + model.objective_function = SQF{T}( + copy(f.quadratic_terms), + copy(f.affine_terms), + f.constant + ) return nothing end -function MOI.get(::QUBOModel{T}, ::MOI.ObjectiveFunctionType) where {T} - return MOI.ScalarQuadraticFunction{T} -end +MOI.get(::QUBOModel{T}, ::MOI.ObjectiveFunctionType) where {T} = SQF{T} function MOI.get(model::QUBOModel, ::MOI.ListOfConstraintTypesPresent) if MOI.is_empty(model) @@ -128,28 +105,35 @@ function MOI.get( model::QUBOModel, ::MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.ZeroOne}, ) - v = MOI.get(model, MOI.ListOfVariableIndices()) - - return [MOI.ConstraintIndex{VI,MOI.ZeroOne}(vi.value) for vi in v] + return [MOI.ConstraintIndex{VI,MOI.ZeroOne}(vi.value) for vi in model.variables] end -function MOI.get( - model::QUBOModel, - ::MOI.ListOfVariableIndices -) - return QUBOTools.variables(model) +function MOI.get(model::QUBOModel, ::MOI.ListOfVariableIndices) + return model.variables end -MOI.get( +function MOI.get( ::QUBOModel, ::MOI.ConstraintFunction, ci::MOI.ConstraintIndex{VI,MOI.ZeroOne}, -) = VI(ci.value) +) + return VI(ci.value) +end -MOI.get( - ::QUBOModel, - ::MOI.ConstraintSet, - ::MOI.ConstraintIndex{VI,MOI.ZeroOne}, -) = MOI.ZeroOne() +function MOI.get(::QUBOModel, ::MOI.ConstraintSet, ::MOI.ConstraintIndex{VI,MOI.ZeroOne}) + return MOI.ZeroOne() +end + +function MOI.get(::QUBOModel, ::MOI.VariableName, vi::VI) + return "x[$(vi.value)]" +end -MOI.get(::QUBOModel, ::MOI.VariableName, vi::VI) = "x[$(vi.value)]" \ No newline at end of file +function MOI.get(::QUBOModel{T}, ::MOI.VariablePrimalStart, ::VI) where {T} + return nothing +end + +MOI.supports(::QUBOModel, ::MOI.VariablePrimalStart, ::MOI.VariableIndex) = true + +function MOI.get(model::QUBOModel, ::MOI.NumberOfVariables) + return length(model.variables) +end diff --git a/src/model/virtual.jl b/src/model/virtual.jl index 165a4110..3d938512 100644 --- a/src/model/virtual.jl +++ b/src/model/virtual.jl @@ -1,13 +1,8 @@ -@doc raw""" - abstract type AbstractVirtualModel{T} <: MOI.AbstractOptimizer end -""" -abstract type AbstractVirtualModel{T} <: MOI.AbstractOptimizer end - -# -*- Virtual Variable Encoding -*- +# Virtual Variable Encoding abstract type Encoding end @doc raw""" - encode!(model::AbstractVirtualModel{T}, v::VirtualVariable{T}) where {T} + encode!(model::VirtualModel{T}, v::VirtualVariable{T}) where {T} Maps newly created virtual variable `v` within the virtual model structure. It follows these steps: @@ -27,111 +22,214 @@ Maps newly created virtual variable `v` within the virtual model structure. It f # References: * [1] Chancellor, N. (2019). Domain wall encoding of discrete variables for quantum annealing and QAOA. _Quantum Science and Technology_, _4_(4), 045004. [{doi}](https://doi.org/10.1088/2058-9565/ab33c2) """ -struct VirtualVariable{T,E<:Encoding} +struct VirtualVariable{T} + e::Encoding x::Union{VI,Nothing} # Source variable (if there is one) y::Vector{VI} # Target variables ξ::PBO.PBF{VI,T} # Expansion function h::Union{PBO.PBF{VI,T},Nothing} # Penalty function (i.e. ‖gᵢ(x)‖ₛ for g(i) ∈ S) - function VirtualVariable{T,E}( + function VirtualVariable{T}( + e::Encoding, x::Union{VI,Nothing}, y::Vector{VI}, ξ::PBO.PBF{VI,T}, h::Union{PBO.PBF{VI,T},Nothing}, - ) where {T,E<:Encoding} - return new{T,E}(x, y, ξ, h) + ) where {T} + return new{T}(e, x, y, ξ, h) end end -# -*- Variable Information -*- +const VV{T} = VirtualVariable{T} + +encoding(v::VirtualVariable) = v.e source(v::VirtualVariable) = v.x target(v::VirtualVariable) = v.y is_aux(v::VirtualVariable) = isnothing(source(v)) expansion(v::VirtualVariable) = v.ξ penaltyfn(v::VirtualVariable) = v.h -# ~*~ Alias ~*~ -const VV{T,E} = VirtualVariable{T,E} +@doc raw""" + VirtualModel{T}(optimizer::Union{Nothing, Type{<:MOI.AbstractOptimizer}} = nothing) where {T} + +This Virtual Model links the final QUBO formulation to the original one, allowing variable value retrieving and other features. +""" +struct VirtualModel{T} <: MOI.AbstractOptimizer + # Underlying Optimizer # + optimizer::Union{MOI.AbstractOptimizer,Nothing} + + # MathOptInterface Bridges # + bridge_model::MOIB.LazyBridgeOptimizer{PreQUBOModel{T}} -function encode!(model::AbstractVirtualModel{T}, v::VV{T}) where {T} + # Virtual Model Interface # + source_model::PreQUBOModel{T} + target_model::QUBOModel{T} + variables::Vector{VV{T}} + source::Dict{VI,VV{T}} + target::Dict{VI,VV{T}} + + # PBO/PBF IR # + f::PBO.PBF{VI,T} # Objective Function + g::Dict{CI,PBO.PBF{VI,T}} # Constraint Functions + h::Dict{VI,PBO.PBF{VI,T}} # Variable Functions + ρ::Dict{CI,T} # Constraint Penalties + θ::Dict{VI,T} # Variable Penalties + H::PBO.PBF{VI,T} # Final Hamiltonian + + # Settings + compiler_settings::Dict{Symbol,Any} + variable_settings::Dict{Symbol,Dict{VI,Any}} + constraint_settings::Dict{Symbol,Dict{CI,Any}} + + function VirtualModel{T}( + constructor::Union{Type{O},Function}; + kws..., + ) where {T,O<:MOI.AbstractOptimizer} + optimizer = constructor() + + return VirtualModel{T}(optimizer; kws...) + end + + function VirtualModel{T}( + optimizer::Union{O,Nothing} = nothing; + kws..., + ) where {T,O<:MOI.AbstractOptimizer} + source_model = PreQUBOModel{T}() + target_model = QUBOModel{T}() + bridge_model = MOIB.full_bridge_optimizer(source_model, T) + + new{T}( + # Underlying Optimizer # + optimizer, + + # MathOptInterface Bridges # + bridge_model, + + # Virtual Model Interface + source_model, + target_model, + Vector{VV{T}}(), + Dict{VI,VV{T}}(), + Dict{VI,VV{T}}(), + + # PBO/PBF IR + PBO.PBF{VI,T}(), # Objective Function + Dict{CI,PBO.PBF{VI,T}}(), # Constraint Functions + Dict{VI,PBO.PBF{VI,T}}(), # Variable Functions + Dict{CI,T}(), # Constraint Penalties + Dict{VI,T}(), # Variable Penalties + PBO.PBF{VI,T}(), # Final Hamiltonian + + # Settings + Dict{Symbol,Any}(), + Dict{Symbol,Dict{VI,Any}}(), + Dict{Symbol,Dict{CI,Any}}(), + ) + end + +end + +VirtualModel(args...; kws...) = VirtualModel{Float64}(args...; kws...) + +function encode!(model::VirtualModel{T}, v::VV{T}) where {T} if !is_aux(v) let x = source(v) - MOI.set(model, Source(), x, v) + model.source[x] = v end end for y in target(v) - MOI.add_constraint(MOI.get(model, TargetModel()), y, MOI.ZeroOne()) - MOI.set(model, Target(), y, v) + MOI.add_constraint(model.target_model, y, MOI.ZeroOne()) + model.target[y] = v end # Add variable to collection - push!(MOI.get(model, Variables()), v) + push!(model.variables, v) return v end -abstract type LinearEncoding <: Encoding end +@doc raw""" + LinearEncoding + +Every linear encoding ``\xi`` is of the form +```math +\xi(\mathbf{y}) = \alpha + \sum_{i = 1}^{n} \gamma_{i} y_{i} +``` + +""" abstract type LinearEncoding <: Encoding end -function VirtualVariable{T,E}( +function VirtualVariable{T}( + e::LinearEncoding, x::Union{VI,Nothing}, y::Vector{VI}, γ::Vector{T}, α::T = zero(T), -) where {T,E<:LinearEncoding} +) where {T} @assert (n = length(y)) == length(γ) ξ = α + PBO.PBF{VI,T}(y[i] => γ[i] for i = 1:n) - return VirtualVariable{T,E}(x, y, ξ, nothing) + return VirtualVariable{T}(e, x, y, ξ, nothing) end function encode!( - ::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::LinearEncoding, x::Union{VI,Nothing}, γ::Vector{T}, α::T = zero(T), -) where {T,E<:LinearEncoding} +) where {T} n = length(γ) - y = MOI.add_variables(MOI.get(model, TargetModel()), n) - v = VirtualVariable{T,E}(x, y, γ, α) + y = MOI.add_variables(model.target_model, n) + v = VirtualVariable{T}(e, x, y, γ, α) return encode!(model, v) end @doc raw""" + Mirror + +Mirrors binary variable ``x \in \mathbb{B}`` with a twin variable ``y \in \mathbb{B}``. """ struct Mirror <: LinearEncoding end -function encode!(e::Mirror, model::AbstractVirtualModel{T}, x::Union{VI,Nothing}) where {T} - return encode!(e, model, x, ones(T, 1)) +function encode!(model::VirtualModel{T}, e::Mirror, x::Union{VI,Nothing}) where {T} + return encode!(model, e, x, ones(T, 1)) end @doc raw""" + Linear """ struct Linear <: LinearEncoding end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::Linear, x::Union{VI,Nothing}, Γ::Function, n::Integer, -) where {T,E<:Linear} +) where {T} γ = T[Γ(i) for i = 1:n] - return encode!(e, model, x, γ, zero(T)) + return encode!(model, e, x, γ, zero(T)) end @doc raw""" + Unary() + +Let ``x \in [a, b] \subset \mathbb{Z}, n = b - a, \mathbf{y} \in \mathbb{B}^{n}``. + +```math +x \approx \xi(\mathbf{y}) = a + \sum_{j = 1}^{b - a} y_{j} +``` """ struct Unary <: LinearEncoding end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::Unary, x::Union{VI,Nothing}, a::T, b::T, -) where {T,E<:Unary} +) where {T} α, β = if a < b ceil(a), floor(b) else @@ -142,37 +240,39 @@ function encode!( M = trunc(Int, β - α) γ = ones(T, M) - return encode!(e, model, x, γ, α) + return encode!(model, e, x, γ, α) end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::Unary, x::Union{VI,Nothing}, a::T, b::T, n::Integer, -) where {T,E<:Unary} +) where {T} Γ = (b - a) / n γ = Γ * ones(T, n) - return encode!(e, model, x, γ, a) + return encode!(model, e, x, γ, a) end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::Unary, x::Union{VI,Nothing}, a::T, b::T, τ::T, -) where {T,E<:Unary} +) where {T} n = ceil(Int, (1 + abs(b - a) / 4τ)) - return encode!(e, model, x, a, b, n) + return encode!(model, e, x, a, b, n) end @doc raw""" + Binary + Binary Expansion within the closed interval ``[\alpha, \beta]``. For a given variable ``x \in [\alpha, \beta]`` we approximate it by @@ -185,12 +285,12 @@ where ``n`` is the number of bits and ``y_i \in \mathbb{B}``. """ struct Binary <: LinearEncoding end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::Binary, x::Union{VI,Nothing}, a::T, b::T, -) where {T,E<:Binary} +) where {T} α, β = if a < b ceil(a), floor(b) else @@ -207,46 +307,47 @@ function encode!( T[[2^i for i = 0:N-2]; [M - 2^(N - 1) + 1]] end - return encode!(e, model, x, γ, α) + return encode!(model, e, x, γ, α) end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::Binary, x::Union{VI,Nothing}, a::T, b::T, n::Integer, -) where {T,E<:Binary} +) where {T} Γ = (b - a) / (2^n - 1) γ = Γ * 2 .^ collect(T, 0:n-1) - return encode!(e, model, x, γ, a) + return encode!(model, e, x, γ, a) end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::Binary, x::Union{VI,Nothing}, a::T, b::T, τ::T, -) where {T,E<:Binary} +) where {T} n = ceil(Int, log2(1 + abs(b - a) / 4τ)) - return encode!(e, model, x, a, b, n) + return encode!(model, e, x, a, b, n) end @doc raw""" + Arithmetic """ struct Arithmetic <: LinearEncoding end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::Arithmetic, x::Union{VI,Nothing}, a::T, b::T, -) where {T,E<:Arithmetic} +) where {T} α, β = if a < b ceil(a), floor(b) else @@ -259,60 +360,77 @@ function encode!( γ = T[[i for i = 1:N-1]; [M - N * (N - 1) / 2]] - return encode!(e, model, x, γ, α) + return encode!(model, e, x, γ, α) end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::Arithmetic, x::Union{VI,Nothing}, a::T, b::T, n::Integer, -) where {T,E<:Arithmetic} +) where {T} Γ = 2 * (b - a) / (n * (n + 1)) γ = Γ * collect(1:n) - return encode!(e, model, x, γ, a) + return encode!(model, e, x, γ, a) end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::Arithmetic, x::Union{VI,Nothing}, a::T, b::T, τ::T, -) where {T,E<:Arithmetic} +) where {T} n = ceil(Int, (1 + sqrt(3 + (b - a) / 2τ)) / 2) - return encode!(e, model, x, a, b, n) + return encode!(model, e, x, a, b, n) end @doc raw""" + OneHot() + +The one-hot encoding is a linear technique used to represent a variable +``x \in \left\lbrace{\gamma_{j}}}\right\rbrace_{j \in \left[n\right]``. + +The encoding function is combined with a constraint assuring that only +one and exactly one of the expansion's variables ``y_{j}`` is activated +at a time. + +```math +\begin{array}{rl} +x = \xi(\mathbf{y}) = & \sum_{j = 1}^{n} \gamma_{j} y_{j} \\ + \mathrm{s.t.} & \sum_{j = 1}^{n} y_{j} = 1 +\end{array} +``` + """ struct OneHot <: LinearEncoding end -function VirtualVariable{T,E}( +function VirtualVariable{T}( + e::OneHot, x::Union{VI,Nothing}, y::Vector{VI}, γ::Vector{T}, α::T = zero(T), -) where {T,E<:OneHot} +) where {T} @assert (n = length(y)) == length(γ) ξ = α + PBO.PBF{VI,T}(y[i] => γ[i] for i = 1:n) h = (one(T) - PBO.PBF{VI,T}(y))^2 - return VirtualVariable{T,E}(x, y, ξ, h) + return VirtualVariable{T}(e, x, y, ξ, h) end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::OneHot, x::Union{VI,Nothing}, a::T, b::T, -) where {T,E<:OneHot} +) where {T} α, β = if a < b ceil(a), floor(b) else @@ -322,75 +440,89 @@ function encode!( # assumes: β - α > 0 γ = collect(T, α:β) - return encode!(e, model, x, γ) + return encode!(model, e, x, γ) end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::OneHot, x::Union{VI,Nothing}, a::T, b::T, n::Integer, -) where {T,E<:OneHot} +) where {T} Γ = (b - a) / (n - 1) γ = a .+ Γ * collect(T, 0:n-1) - return encode!(e, model, x, γ) + return encode!(model, e, x, γ) end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::OneHot, x::Union{VI,Nothing}, a::T, b::T, τ::T, -) where {T,E<:OneHot} +) where {T} n = ceil(Int, (1 + abs(b - a) / 4τ)) - return encode!(e, model, x, a, b, n) + return encode!(model, e, x, a, b, n) end abstract type SequentialEncoding <: Encoding end function encode!( - ::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::SequentialEncoding, x::Union{VI,Nothing}, γ::Vector{T}, α::T = zero(T), -) where {T,E<:SequentialEncoding} +) where {T} n = length(γ) - y = MOI.add_variables(MOI.get(model, TargetModel()), n - 1) - v = VirtualVariable{T,E}(x, y, γ, α) + y = MOI.add_variables(model.target_model, n - 1) + v = VirtualVariable{T}(e, x, y, γ, α) return encode!(model, v) end -struct DomainWall <: SequentialEncoding end +@doc raw""" + DomainWall() + +The Domain Wall[^Chancellor2019] encoding method is a sequential approach that requires only +``n - 1`` bits to represent ``n`` distinct values. + +!!! table "Encoding Analysis" + | | bits | linear | quadratic | ``\Delta`` | + | :-: | :--: | :----: | :-------: | :--------: | + | Domain Wall | ``n - 1`` | ``n`` | | ``O(n)`` | + +[^Chancellor2019]: + Nicholas Chancellor, **Domain wall encoding of discrete variables for quantum annealing and QAOA**, *Quantum Science Technology 4*, 2019. +""" struct DomainWall <: SequentialEncoding end -function VirtualVariable{T,E}( +function VirtualVariable{T}( + e::DomainWall, x::Union{VI,Nothing}, y::Vector{VI}, γ::Vector{T}, α::T = zero(T), -) where {T,E<:DomainWall} +) where {T} @assert (n = length(y)) == length(γ) - 1 ξ = α + PBO.PBF{VI,T}(y[i] => (γ[i] - γ[i+1]) for i = 1:n) h = 2 * (PBO.PBF{VI,T}(y[2:n]) - PBO.PBF{VI,T}([Set{VI}([y[i], y[i-1]]) for i = 2:n])) - return VirtualVariable{T,E}(x, y, ξ, h) + return VirtualVariable{T}(e, x, y, ξ, h) end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::DomainWall, x::Union{VI,Nothing}, a::T, b::T, -) where {T,E<:DomainWall} +) where {T} α, β = if a < b ceil(a), floor(b) else @@ -401,252 +533,164 @@ function encode!( M = trunc(Int, β - α) γ = α .+ collect(T, 0:M) - return encode!(e, model, x, γ) + return encode!(model, e, x, γ) end function encode!( - e::E, - model::AbstractVirtualModel{T}, + model::VirtualModel{T}, + e::DomainWall, x::Union{VI,Nothing}, a::T, b::T, n::Integer, -) where {T,E<:DomainWall} +) where {T} Γ = (b - a) / (n - 1) γ = a .+ Γ * collect(T, 0:n-1) - return encode!(e, model, x, γ) -end - -mutable struct VirtualQUBOModelSettings{T} - atol::Dict{Union{VI,Nothing},T} - bits::Dict{Union{VI,Nothing},Int} - encoding::Dict{Union{CI,VI,Nothing},Encoding} - - function VirtualQUBOModelSettings{T}(; - atol::T = 1E-2, - bits::Integer = 3, - encoding::Encoding = Binary(), - ) where {T} - return new{T}( - Dict{Union{VI,Nothing},T}(nothing => atol), - Dict{Union{VI,Nothing},Int}(nothing => bits), - Dict{Union{VI,Nothing},Encoding}(nothing => encoding), - ) - end + return encode!(model, e, x, γ) end @doc raw""" - VirtualQUBOModel{T}(optimizer::Union{Nothing, Type{<:MOI.AbstractOptimizer}} = nothing) where {T} - -This QUBO Virtual Model links the final QUBO formulation to the original one, allowing variable value retrieving and other features. -""" -struct VirtualQUBOModel{T} <: AbstractVirtualModel{T} - # -*- Underlying Optimizer -*- # - optimizer::Union{MOI.AbstractOptimizer,Nothing} - - # -*- MathOptInterface Bridges -*- # - bridge_model::MOIB.LazyBridgeOptimizer{PreQUBOModel{T}} - - # -*- Virtual Model Interface -*- # - source_model::PreQUBOModel{T} - target_model::QUBOModel{T} - variables::Vector{VV{T}} - source::Dict{VI,VV{T}} - target::Dict{VI,VV{T}} + Bounded{Binary,T}(μ::T) where {T} - # -*- PBO/PBF IR -*- # - f::PBO.PBF{VI,T} # Objective Function - g::Dict{CI,PBO.PBF{VI,T}} # Constraint Functions - h::Dict{VI,PBO.PBF{VI,T}} # Variable Functions - ρ::Dict{CI,T} # Constraint Penalties - θ::Dict{VI,T} # Variable Penalties - H::PBO.PBF{VI,T} # Final Hamiltonian - - # -*- Settings -*- - settings::VirtualQUBOModelSettings{T} +Let ``x \in [a, b] \subset \mathbb{Z}`` and ``n = b - a``. - function VirtualQUBOModel{T}( - constructor::Union{Type{O},Function}; - kws..., - ) where {T,O<:MOI.AbstractOptimizer} - optimizer = constructor() - - return VirtualQUBOModel{T}(optimizer; kws...) - end - - function VirtualQUBOModel{T}( - optimizer::Union{O,Nothing} = nothing; - kws..., - ) where {T,O<:MOI.AbstractOptimizer} - source_model = PreQUBOModel{T}() - target_model = QUBOModel{T}() - bridge_model = MOIB.full_bridge_optimizer(source_model, T) - - new{T}( - # -*- Underlying Optimizer -*- # - optimizer, - - # -*- MathOptInterface Bridges -*- # - bridge_model, - - # -*- Virtual Model Interface -*- - source_model, - target_model, - Vector{VV{T}}(), - Dict{VI,VV{T}}(), - Dict{VI,VV{T}}(), - - # -*- PBO/PBF IR -*- - PBO.PBF{VI,T}(), # Objective Function - Dict{CI,PBO.PBF{VI,T}}(), # Constraint Functions - Dict{VI,PBO.PBF{VI,T}}(), # Variable Functions - Dict{CI,T}(), # Constraint Penalties - Dict{VI,T}(), # Variable Penalties - PBO.PBF{VI,T}(), # Final Hamiltonian +First, - # -*- Settings -*- - VirtualQUBOModelSettings{T}(; kws...), - ) - end - - VirtualQUBOModel(args...; kws...) = VirtualQUBOModel{Float64}(args...; kws...) -end - -QUBOTools.backend(model::VirtualQUBOModel) = QUBOTools.backend(model.target_model) - -struct Source <: MOI.AbstractModelAttribute end -function MOI.get(::AbstractVirtualModel, ::Source) end -function MOI.get(::AbstractVirtualModel, ::Source, ::VI) end -function MOI.set(::AbstractVirtualModel, ::Source, ::VI, ::VV) end +```math +\begin{align*} + 2^{k - 1} &\le \mu \\ +\implies k &= \left\lfloor\log_{2} \mu + 1 \right\rfloor +\end{align*} +``` -struct Target <: MOI.AbstractModelAttribute end -function MOI.get(::AbstractVirtualModel, ::Target) end -function MOI.get(::AbstractVirtualModel, ::Target, ::VI) end -function MOI.set(::AbstractVirtualModel, ::Target, ::VI, ::VV) end +Since -struct Variables <: MOI.AbstractModelAttribute end -function MOI.get(::AbstractVirtualModel, ::Variables) end +```math +\sum_{j = 1}^{k} 2^{j - 1} = \sum_{j = 0}^{k - 1} 2^{j} = 2^{k} - 1 +``` -struct SourceModel <: MOI.AbstractModelAttribute end -function MOI.get(::AbstractVirtualModel, ::SourceModel) end +then, for ``r \in \mathbb{N}`` -struct TargetModel <: MOI.AbstractModelAttribute end -function MOI.get(::AbstractVirtualModel, ::TargetModel) end +```math +n = 2^{k} - 1 + r \times \mu + \epsilon \implies r = \left\lfloor \frac{n - 2^{k} + 1}{\mu} \right\rfloor +``` -function MOI.is_empty(model::AbstractVirtualModel) - return all([ - MOI.is_empty(MOI.get(model, SourceModel())), - MOI.is_empty(MOI.get(model, TargetModel())), - ]) -end +and -function MOI.empty!(model::AbstractVirtualModel) - MOI.empty!(MOI.get(model, SourceModel())) - MOI.empty!(MOI.get(model, TargetModel())) - empty!(MOI.get(model, Variables())) +```math +\epsilon = n - 2^{k} + 1 - r \times \mu +``` - return nothing -end +\gamma_{j} = \left\lbrace\begin{array}{cl} + 2^{j} & \text{if } 1 \le j \le k \\ + \mu & \text{if } k < j < r + k \\ + n - 2^k + 1 - r \times \mu & \text{otherwise} +\end{array}\right. +``` -function MOI.get( - model::AbstractVirtualModel, - attr::Union{ - MOI.ListOfConstraintAttributesSet, - MOI.ListOfConstraintIndices, - MOI.ListOfConstraintTypesPresent, - MOI.ListOfModelAttributesSet, - MOI.ListOfVariableAttributesSet, - MOI.ListOfVariableIndices, - MOI.NumberOfConstraints, - MOI.NumberOfVariables, - MOI.Name, - MOI.ObjectiveFunction, - MOI.ObjectiveFunctionType, - MOI.ObjectiveSense, - }, -) - - return MOI.get(MOI.get(model, SourceModel()), attr) -end + Bounded{Unary,T}(μ::T) where {T} -function MOI.get( - model::AbstractVirtualModel, - attr::Union{MOI.ConstraintFunction,MOI.ConstraintSet}, - ci::MOI.ConstraintIndex, -) - return MOI.get(MOI.get(model, SourceModel()), attr, ci) -end +""" struct Bounded{E<:LinearEncoding,T} <: LinearEncoding + μ::T -function MOI.set( - model::AbstractVirtualModel, - attr::Union{MOI.ObjectiveFunction,MOI.ObjectiveSense}, - value::Any, -) + function Bounded{E,T}(μ::T) where {E,T} + @assert !iszero(μ) - return MOI.get(MOI.get(model, SourceModel()), attr, value) + return new{E,T}(μ) + end end -function MOI.get(model::AbstractVirtualModel, attr::MOI.VariableName, x::VI) - return MOI.get(MOI.get(model, SourceModel()), attr, x) +function Bounded{E}(μ::T) where {E,T} + return Bounded{E,T}(μ) end -MOI.get(model::VirtualQUBOModel, ::Source) = model.source -MOI.get(model::VirtualQUBOModel, ::Source, x::VI) = model.source[x] +function encode!( + model::VirtualModel{T}, + e::Bounded{Binary,T}, + x::Union{VI,Nothing}, + a::T, + b::T, +) where {T} + if a < b + a = ceil(a) + b = floor(b) + else + a = ceil(b) + b = floor(a) + end -function MOI.set(model::VirtualQUBOModel{T}, ::Source, x::VI, v::VV{T}) where {T} - model.source[x] = v -end + n = round(Int, b - a) + k = floor(Int, log2(e.μ) + 1) + m = 2^k - 1 + r = floor(Int, (n - m) / e.μ) + ϵ = n - m - r * e.μ -MOI.get(model::VirtualQUBOModel, ::Target) = model.target -MOI.get(model::VirtualQUBOModel, ::Target, y::VI) = model.target[y] + if iszero(ϵ) + γ = T[[2^(j - 1) for j = 1:k]; [e.μ for _ = 1:r]] + else + γ = T[[2^(j - 1) for j = 1:k]; [e.μ for _ = 1:r]; [ϵ]] + end -function MOI.set(model::VirtualQUBOModel{T}, ::Target, y::VI, v::VV{T}) where {T} - model.target[y] = v + return encode!(model, e, x, γ, a) end -MOI.get(model::VirtualQUBOModel, ::Variables) = model.variables -MOI.get(model::VirtualQUBOModel, ::SourceModel) = model.source_model -MOI.get(model::VirtualQUBOModel, ::TargetModel) = model.target_model - -function Base.show(io::IO, model::AbstractVirtualModel) - print( - io, - """ - Virtual Model - with source: - $(MOI.get(model, SourceModel())) - with target: - $(MOI.get(model, TargetModel())) - """, - ) -end +function encode!( + model::VirtualModel{T}, + e::Bounded{Unary,T}, + x::Union{VI,Nothing}, + a::T, + b::T, +) where {T} + if a < b + a = ceil(a) + b = floor(b) + else + a = ceil(b) + b = floor(a) + end -function MOI.add_variable(model::VirtualQUBOModel) - source_model = MOI.get(model, SourceModel()) + n = round(Int, b - a) + k = ceil(Int, e.μ - 1) + r = floor(Int, (n - k) / e.μ) + ϵ = n - k + - r * e.μ - return MOI.add_variable(source_model) -end + if iszero(ϵ) + γ = T[ones(T, k); [e.μ for _ = 1:r]] + else + γ = T[ones(T, k); [e.μ for _ = 1:r]; [ϵ]] + end -function MOI.add_constraint( - model::VirtualQUBOModel, - f::MOI.AbstractFunction, - s::MOI.AbstractSet, -) - source_model = MOI.get(model, SourceModel()) - - return MOI.add_constraint(source_model, f, s) + return encode!(model, e, x, γ, a) end -function MOI.set(model::VirtualQUBOModel, os::MOI.ObjectiveSense, s::MOI.OptimizationSense) - source_model = MOI.get(model, SourceModel()) +function encode!( + model::VirtualModel{T}, + e::Bounded{Arithmetic,T}, + x::Union{VI,Nothing}, + a::T, + b::T, +) where {T} + if a < b + a = ceil(a) + else + b = floor(b) + a = ceil(b) + b = floor(a) + end - MOI.set(source_model, os, s) -end + n = round(Int, b - a) + k = floor(Int, e.μ) + m = (k * (k + 1)) ÷ 2 + r = floor(Int, (n - m) / e.μ) + ϵ = n - m + - r * e.μ -function MOI.set(model::VirtualQUBOModel, of::MOI.ObjectiveFunction, f::MOI.AbstractFunction) - source_model = MOI.get(model, SourceModel()) + if iszero(ϵ) + γ = T[collect(T,1:k); [e.μ for _ = 1:r]] + else + γ = T[collect(T,1:k); [e.μ for _ = 1:r]; [ϵ]] + end - MOI.set(source_model, of, f) + return encode!(model, e, x, γ, a) end diff --git a/src/model/wrapper.jl b/src/model/wrapper.jl index 6b570dad..bab0f250 100644 --- a/src/model/wrapper.jl +++ b/src/model/wrapper.jl @@ -1,188 +1,106 @@ -function MOI.empty!(model::VirtualQUBOModel) - # -*- Models -*- - MOI.empty!(MOI.get(model, SourceModel())) - MOI.empty!(MOI.get(model, TargetModel())) - - # -*- Virtual Variables -*- - empty!(MOI.get(model, Variables())) - empty!(MOI.get(model, Source())) - empty!(MOI.get(model, Target())) - - # -*- Underlying Optimizer -*- - if !isnothing(model.optimizer) - MOI.empty!(model.optimizer) - end - - # -*- PBF/IR -*- - empty!(model.f) - empty!(model.g) - empty!(model.h) - empty!(model.ρ) - empty!(model.θ) - - return nothing -end - # Notes on the optimize! interface # After `JuMP.optimize!(model)` there are a few layers before reaching -# 1. `MOI.optimize!(::VirtualQUBOModel, ::MOI.ModelLike)` +# 1. `MOI.optimize!(::VirtualModel, ::MOI.ModelLike)` # Then, -# 2. `MOI.copy_to(::VirtualQUBOModel, ::MOI.ModelLike)` -# 3. `MOI.optimize!(::VirtualQUBOModel)` +# 2. `MOI.copy_to(::VirtualModel, ::MOI.ModelLike)` +# 3. `MOI.optimize!(::VirtualModel)` # is called. -function MOI.optimize!(model::VirtualQUBOModel) - source_model = MOI.get(model, SourceModel()) - target_model = MOI.get(model, TargetModel()) - index_map = MOIU.identity_index_map(source_model) +function MOI.optimize!(model::VirtualModel) + index_map = MOIU.identity_index_map(model.source_model) - # -*- JuMP to QUBO Compilation -*- # + # De facto JuMP to QUBO Compilation ToQUBO.toqubo!(model) if !isnothing(model.optimizer) - MOI.optimize!(model.optimizer, target_model) + MOI.optimize!(model.optimizer, model.target_model) end return (index_map, false) end -function MOI.copy_to(model::VirtualQUBOModel{T}, source::MOI.ModelLike) where {T} +function MOI.copy_to(model::VirtualModel{T}, source::MOI.ModelLike) where {T} if !MOI.is_empty(model) error("QUBO Model is not empty") end - # -*- Copy to PreQUBOModel + Add Bridges -*- # - source_model = MOI.get(model, SourceModel()) - bridge_model = MOIB.full_bridge_optimizer(source_model, T) + # Copy to PreQUBOModel + Add Bridges + bridge_model = MOIB.full_bridge_optimizer(model.source_model, T) - # -*- Copy to source using bridges - *- # + # Copy to source using bridges return MOI.copy_to(bridge_model, source) # index_map end -# -*- :: Objective Function Support :: -*- # +# Objective Function Support MOI.supports( - ::VirtualQUBOModel{T}, + ::VirtualModel{T}, ::MOI.ObjectiveFunction{<:Union{VI,SAF{T},SQF{T}}}, ) where {T} = true -# -*- :: Constraint Support :: -*- # +# Constraint Support MOI.supports_constraint( - ::VirtualQUBOModel{T}, + ::VirtualModel{T}, ::Type{VI}, - ::Type{<:Union{ - MOI.ZeroOne, - MOI.Integer, - MOI.Interval{T}, - MOI.LessThan{T}, - MOI.GreaterThan{T}, - }}, + ::Type{ + <:Union{MOI.ZeroOne,MOI.Integer,MOI.Interval{T},MOI.LessThan{T},MOI.GreaterThan{T}}, + }, ) where {T} = true MOI.supports_constraint( - ::VirtualQUBOModel{T}, + ::VirtualModel{T}, ::Type{<:Union{SAF{T},SQF{T}}}, ::Type{<:Union{MOI.EqualTo{T},MOI.LessThan{T}}}, ) where {T} = true MOI.supports_constraint( - ::VirtualQUBOModel{T}, + ::VirtualModel{T}, ::Type{<:MOI.VectorOfVariables}, ::Type{<:MOI.SOS1}, ) where {T} = true MOI.supports_add_constrained_variable( - ::VirtualQUBOModel{T}, - ::Type{<:Union{ - MOI.ZeroOne, - MOI.Integer, - MOI.Interval{T}, - MOI.LessThan{T}, - MOI.GreaterThan{T}, - }}, + ::VirtualModel{T}, + ::Type{ + <:Union{MOI.ZeroOne,MOI.Integer,MOI.Interval{T},MOI.LessThan{T},MOI.GreaterThan{T}}, + }, ) where {T} = true -function MOI.get( - model::VirtualQUBOModel, - attr::Union{ - MOI.SolveTimeSec, - MOI.PrimalStatus, - MOI.DualStatus, - MOI.TerminationStatus, - MOI.RawStatusString, - } -) - if !isnothing(model.optimizer) - MOI.get(model.optimizer, attr) - else - nothing - end -end -function MOI.set( - model::VirtualQUBOModel, - attr::Union{ - MOI.SolveTimeSec, - MOI.PrimalStatus, - MOI.DualStatus, - MOI.TerminationStatus, - MOI.RawStatusString, - }, - value::Any, -) - return MOI.set(model.optimizer, attr, value) -end +PBO.showvar(x::VI) = PBO.showvar(x.value) -MOI.supports(::VirtualQUBOModel, ::MOI.SolveTimeSec) = true +PBO.varlt(x::VI, y::VI) = PBO.varlt(x.value, y.value) -function MOI.get(model::VirtualQUBOModel, rc::MOI.ResultCount) - if isnothing(model.optimizer) - return 0 - else - return MOI.get(model.optimizer, rc) - end -end +function PBO.varlt(x::Set{V}, y::Set{V}) where {V} + if length(x) == length(y) + xv = sort!(collect(x); lt = PBO.varlt) + yv = sort!(collect(y); lt = PBO.varlt) -MOI.supports(::VirtualQUBOModel, ::MOI.ResultCount) = true + for (xi, yi) in zip(xv, yv) + if xi == yi + continue + else + return PBO.varlt(xi, yi) + end + end -function MOI.get(model::VirtualQUBOModel{T}, ov::MOI.ObjectiveValue) where {T} - if isnothing(model.optimizer) - return zero(T) + return false else - return MOI.get(model.optimizer, ov) + return length(x) < length(y) end end -function MOI.get(model::VirtualQUBOModel{T}, vp::MOI.VariablePrimal, x::VI) where {T} - if isnothing(model.optimizer) - return zero(T) - else - v = MOI.get(model, Source(), x) - s = zero(T) +const Optimizer{T} = VirtualModel{T} - for (ω, c) in expansion(v) - for y in ω - c *= MOI.get(model.optimizer, vp, y) - end - - s += c - end - - return s - end -end - -MOI.get(::VirtualQUBOModel, ::MOI.SolverName) = "Virtual QUBO Model" -MOI.get(::VirtualQUBOModel, ::MOI.SolverVersion) = PROJECT_VERSION +# QUBOTools +function qubo(model, type::Type = Dict) + n, L, Q, α, β = MOI.get(model, Attributes.QUBONormalForm()) -function MOI.get(model::VirtualQUBOModel, rs::MOI.RawSolver) - if isnothing(model.optimizer) - return nothing - else - return MOI.get(model.optimizer, rs) - end + return QUBOTools.qubo(type, n, L, Q, α, β) end -PBO.showvar(x::VI) = PBO.showvar(x.value) -PBO.varcmp(x::VI, y::VI) = PBO.varcmp(x.value, y.value) +function ising(model, type::Type = Dict) + n, L̄, Q̄, ᾱ, β̄ = MOI.get(model, Attributes.QUBONormalForm()) + L, Q, α, β = QUBOTools.cast(QUBOTools.𝔹, QUBOTools.𝕊, L̄, Q̄, ᾱ, β̄) -const Optimizer{T} = VirtualQUBOModel{T} + return QUBOTools.ising(type, n, L, Q, α, β) +end diff --git a/test/Project.toml b/test/Project.toml index a1d40796..fed2a44f 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -7,5 +7,5 @@ TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] -Anneal = "0.5" -JuMP = "1" +Anneal = "0.6.1" +JuMP = "1" diff --git a/test/examples/continuous/continuous.jl b/test/examples/continuous/continuous.jl new file mode 100644 index 00000000..c21a4d47 --- /dev/null +++ b/test/examples/continuous/continuous.jl @@ -0,0 +1,7 @@ +include("continuous_1.jl") + +function test_continuous() + @testset "Continuous Variables" verbose = true begin + test_continuous_1() + end +end \ No newline at end of file diff --git a/test/examples/continuous/continuous_1.jl b/test/examples/continuous/continuous_1.jl new file mode 100644 index 00000000..1641de33 --- /dev/null +++ b/test/examples/continuous/continuous_1.jl @@ -0,0 +1,67 @@ +function test_continuous_1() + @testset "9 variables ∈ [0, 1]" begin + # Problem Data + n = 3 + A = [ + -1 2 2 + 2 -1 2 + 2 2 -1 + ] + + # Solution + Q̄ = [ + -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 + ] + + ᾱ = 1 + β̄ = 0 + + x̄ = Set{Matrix{Int}}([ + [0 1 1;1 0 1;1 1 0] + ]) + ȳ = 12 + + # Model + model = Model(() -> ToQUBO.Optimizer(ExactSampler.Optimizer)) + + set_optimizer_attribute(model, TQA.DefaultVariableEncodingATol(), 1E-1) + + @variable(model, 0 <= x[1:n,1:n] <= 1) + @objective(model, Max, sum(A .* x)) + + optimize!(model) + + # Reformulation + Q, α, β = ToQUBO.qubo(model, Matrix) + + @test α ≈ ᾱ + @test β ≈ β̄ + @test Q ≈ Q̄/3 + + # Solutions + x̂ = trunc.(Int, value.(x)) + ŷ = objective_value(model) + + @test x̂ ∈ x̄ + @test ŷ ≈ ȳ + + return nothing + end +end \ No newline at end of file diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 6fc1c52f..b515d813 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -3,6 +3,7 @@ include("linear/linear.jl") include("quadratic/quadratic.jl") include("logical/logical.jl") include("integer/integer.jl") +include("continuous/continuous.jl") function test_examples() @testset "Examples" verbose = true begin @@ -11,5 +12,6 @@ function test_examples() test_quadratic() test_logical() test_integer() + test_continuous() end end \ No newline at end of file diff --git a/test/examples/integer/integer_primes.jl b/test/examples/integer/integer_primes.jl index 35cbb07a..5ff1ade9 100644 --- a/test/examples/integer/integer_primes.jl +++ b/test/examples/integer/integer_primes.jl @@ -1,24 +1,24 @@ function test_integer_primes() @testset "Prime Factoring: 15 = 3 × 5" begin - # ~*~ Problem Data ~*~ # + # Problem Data # R = 15 a = ceil(Int, √R) b = ceil(Int, R ÷ 2) - # ~*~ Solution Data ~*~ # + # Solution Data # ᾱ = 1 β̄ = 49 Q̄ = [ - -24 16 15 15 8 0 20 18 0 20 - 0 -40 60 60 -8 8 -20 0 40 -20 - 0 0 -40 98 -8 -8 -20 -18 -40 0 - 0 0 0 -40 -8 -8 0 -18 -40 -20 - 0 0 0 0 16 0 0 0 0 0 - 0 0 0 0 0 8 0 0 0 0 - 0 0 0 0 0 0 20 0 0 0 - 0 0 0 0 0 0 0 18 0 0 - 0 0 0 0 0 0 0 0 40 0 - 0 0 0 0 0 0 0 0 0 20 + -24 16 15 15 20 20 18 0 8 0 + 0 -40 60 60 -20 -20 0 40 -8 8 + 0 0 -40 98 -20 0 -18 -40 -8 -8 + 0 0 0 -40 0 -20 -18 -40 -8 -8 + 0 0 0 0 20 0 0 0 0 0 + 0 0 0 0 0 20 0 0 0 0 + 0 0 0 0 0 0 18 0 0 0 + 0 0 0 0 0 0 0 40 0 0 + 0 0 0 0 0 0 0 0 16 0 + 0 0 0 0 0 0 0 0 0 8 ] @@ -30,20 +30,22 @@ function test_integer_primes() @variable(model, 2 <= p <= a, Int) @variable(model, a <= q <= b, Int) - @constraint(model, c1, p * q - R == 0) + @constraint(model, c1, p * q == R) + + set_optimizer_attribute(model, TQA.StableQuadratization(), true) optimize!(model) - # :: Reformulation :: # - ρ = MOI.get(unsafe_backend(model), ToQUBO.Penalty(), c1.index) - Q, α, β = ToQUBO.qubo(unsafe_backend(model), Matrix) + # Reformulation + ρ = MOI.get(model, TQA.ConstraintEncodingPenalty(), c1) + Q, α, β = ToQUBO.qubo(model, Matrix) @test ρ ≈ ρ̄ @test α ≈ ᾱ @test β ≈ β̄ - @test Q ≈ Q̄ broken = true + @test Q ≈ Q̄ - # :: Solutions :: # + # Solutions p̂ = trunc(Int, value(p)) q̂ = trunc(Int, value(q)) diff --git a/test/examples/linear/linear1.jl b/test/examples/linear/linear1.jl index e2b00aca..4fb679f8 100644 --- a/test/examples/linear/linear1.jl +++ b/test/examples/linear/linear1.jl @@ -1,6 +1,6 @@ function test_linear1() @testset "3 variables, 1 constraint" begin - # ~*~ Problem Data ~*~ # + # Problem Data n = 3 a = [0.3, 0.5, 1.0] b = 1.6 @@ -9,7 +9,7 @@ function test_linear1() # Penalty Choice ρ̄ = -7 - # ~*~ Solution Data ~*~ # + # Solution Data Q̄ = [ 610 -210 -420 -42 -84 -168 -336 -42 0 947 -700 -70 -140 -280 -560 -70 @@ -27,7 +27,7 @@ function test_linear1() x̄ = Set{Vector{Int}}([[0, 1, 1]]) ȳ = 5 - # ~*~ Model ~*~ # + # Model model = Model(() -> ToQUBO.Optimizer(ExactSampler.Optimizer)) @variable(model, x[1:n], Bin) @@ -36,22 +36,22 @@ function test_linear1() optimize!(model) - # :: Reformulation :: - qubo_model = unsafe_backend(model) - - ρ = MOI.get(qubo_model, ToQUBO.Penalty(), c1.index) - Q, α, β = ToQUBO.qubo(qubo_model, Matrix) + # Reformulation + ρ = MOI.get(model, TQA.ConstraintEncodingPenalty(), c1) + Q, α, β = ToQUBO.qubo(model, Matrix) @test ρ ≈ ρ̄ @test α ≈ ᾱ @test β ≈ β̄ @test Q ≈ Q̄ - # :: Solutions :: + # Solutions x̂ = trunc.(Int, value.(x)) ŷ = objective_value(model) @test x̂ ∈ x̄ @test ŷ ≈ ȳ + + return nothing end end \ No newline at end of file diff --git a/test/examples/linear/linear2.jl b/test/examples/linear/linear2.jl index 486a421e..6d9cabd7 100644 --- a/test/examples/linear/linear2.jl +++ b/test/examples/linear/linear2.jl @@ -1,15 +1,16 @@ function test_linear2() @testset "11 variables, 3 constraints" begin - # ~*~ Problem Data ~*~ # + # Problem Data + m = 3 n = 11 A = [1 0 0 1 1 1 0 1 1 1 1; 0 1 0 1 0 1 1 0 1 1 1; 0 0 1 0 1 0 1 1 1 1 1] b = [1, 1, 1] c = [2, 4, 4, 4, 4, 4, 5, 4, 5, 6, 5] # Penalty Choice - ρ̄ = sum(abs.(c)) + 1 + ρ̄ = fill(sum(abs.(c)) + 1, m) - # ~*~ Solution Data ~*~ # + # Solution Data Q̄ = [ -46 0 0 96 96 96 0 96 96 96 96 0 -44 0 96 0 96 96 0 96 96 96 @@ -34,7 +35,7 @@ function test_linear2() ȳ = 5 - # ~*~ Model ~*~ # + # Model model = Model(() -> ToQUBO.Optimizer(ExactSampler.Optimizer)) @variable(model, x[1:n], Bin) @@ -43,24 +44,22 @@ function test_linear2() optimize!(model) - # :: Reformulation :: - qubo_model = unsafe_backend(model) + # Reformulation + ρ = MOI.get.(model, TQA.ConstraintEncodingPenalty(), k) + Q, α, β = ToQUBO.qubo(model, Matrix) - k = [ki.index for ki in k] - ρ = MOI.get.(qubo_model, ToQUBO.Penalty(), k) - Q, α, β = ToQUBO.PBO.qubo(qubo_model, Matrix) - - @test all(ρ .≈ ρ̄) - + @test ρ ≈ ρ̄ @test α ≈ ᾱ @test β ≈ β̄ @test Q ≈ Q̄ - # :: Solutions :: + # Solutions x̂ = trunc.(Int, value.(x)) ŷ = objective_value(model) @test x̂ ∈ x̄ @test ŷ ≈ ȳ + + return nothing end end \ No newline at end of file diff --git a/test/examples/logical/logical_sos.jl b/test/examples/logical/logical_sos.jl index 571239dd..bc6f4a51 100644 --- a/test/examples/logical/logical_sos.jl +++ b/test/examples/logical/logical_sos.jl @@ -1,6 +1,6 @@ function test_logical_sos1() @testset "SOS1: 3 variables" begin - # ~*~ Problem Data ~*~ # + # Problem Data n = 3 A = [ -1 2 2 @@ -11,7 +11,7 @@ function test_logical_sos1() # Penalty Choice ρ̄ = -16 - # ~*~ Solution Data ~*~ # + # Solution Data Q̄ = [ 15 -28 -28 -32 0 15 -28 -32 @@ -25,7 +25,7 @@ function test_logical_sos1() x̄ = Set{Vector{Int}}([[0, 0, 0]]) ȳ = 0 - # ~*~ Model ~*~ # + # Model model = Model(() -> ToQUBO.Optimizer(ExactSampler.Optimizer)) @variable(model, x[1:n], Bin) @@ -34,22 +34,22 @@ function test_logical_sos1() optimize!(model) - # :: Reformulation :: - qubo_model = unsafe_backend(model) - - ρ = MOI.get(qubo_model, ToQUBO.Penalty(), c1.index) - Q, α, β = ToQUBO.qubo(qubo_model, Matrix) + # Reformulation + ρ = MOI.get(model, TQA.ConstraintEncodingPenalty(), c1) + Q, α, β = ToQUBO.qubo(model, Matrix) @test ρ ≈ ρ̄ @test α ≈ ᾱ @test β ≈ β̄ @test Q ≈ Q̄ - # :: Solutions :: + # Solutions x̂ = trunc.(Int, value.(x)) ŷ = objective_value(model) @test x̂ ∈ x̄ @test ŷ ≈ ȳ + + return nothing end end \ No newline at end of file diff --git a/test/examples/logical/logical_tsp.jl b/test/examples/logical/logical_tsp.jl index 84e0981d..a40f8f5e 100644 --- a/test/examples/logical/logical_tsp.jl +++ b/test/examples/logical/logical_tsp.jl @@ -1,59 +1,73 @@ function test_logical_tsp() @testset "TSP: 16 variables" begin - # ~*~ Problem Data ~*~ # + # Problem Data # n = 4 D = [ - 0 10 10 10 - 10 0 10 10 - 10 10 0 10 - 10 10 10 0 + 0 1 5 4 + 1 0 2 6 + 5 2 0 3 + 4 6 3 0 ] # Penalty Choice - ρ̄ = -16 + ρ̄ = fill(169, 2n) - # ~*~ Solution Data ~*~ # + # Solution Data Q̄ = [ - 15 -28 -28 -32 - 0 15 -28 -32 - 0 0 15 -32 - 0 0 0 16 + -338 676 1 676 5 0 676 5 5 4 + 0 -675 676 681 679 5 681 682 6 6 + 0 0 -338 3 676 2 6 676 6 0 + 0 0 0 -676 681 676 681 10 681 7 + 0 0 0 0 -674 676 4 682 681 6 + 0 0 0 0 0 -338 5 5 676 3 + 0 0 0 0 0 0 -672 682 683 676 + 0 0 0 0 0 0 0 -676 682 676 + 0 0 0 0 0 0 0 0 -673 676 + 0 0 0 0 0 0 0 0 0 -338 ] ᾱ = 1 - β̄ = -16 + β̄ = 1352 + x̄ = Set{Matrix{Int}}([ + [0 0 0 1; 0 0 1 0; 0 1 0 0; 1 0 0 0], + [1 0 0 0; 0 0 0 1; 0 0 1 0; 0 1 0 0], + [0 1 0 0; 1 0 0 0; 0 0 0 1; 0 0 1 0], + [0 0 1 0; 0 1 0 0; 1 0 0 0; 0 0 0 1], + ]) + ȳ = 10 - x̄ = Set{Vector{Int}}([[0, 0, 0]]) - ȳ = 0 - - # ~*~ Model ~*~ # + # Model # model = Model(() -> ToQUBO.Optimizer(ExactSampler.Optimizer)) @variable(model, x[1:n, 1:n], Bin, Symmetric) - @objective(model, Min, sum(D[i, j] * x[i, k] * x[j, (k % n) + 1] for i = 1:n, j = 1:n, k = 1:n)) + @objective( + model, + Min, + sum(D[i, j] * x[i, k] * x[j, (k%n)+1] for i = 1:n, j = 1:n, k = 1:n) + ) @constraint(model, ci[i = 1:n], sum(x[i, :]) == 1) @constraint(model, ck[k = 1:n], sum(x[:, k]) == 1) optimize!(model) - # :: Reformulation :: - qubo_model = unsafe_backend(model) - - ρi = MOI.get.(qubo_model, ToQUBO.Penalty(), collect(map(i->i.index, ci))) - ρk = MOI.get.(qubo_model, ToQUBO.Penalty(), collect(map(i->i.index, ck))) + # Reformulation + ρi = MOI.get.(model, TQA.ConstraintEncodingPenalty(), ci) + ρk = MOI.get.(model, TQA.ConstraintEncodingPenalty(), ck) ρ = [ρi; ρk] - Q, α, β = ToQUBO.qubo(qubo_model, Matrix) + Q, α, β = ToQUBO.qubo(model, Matrix) - @show ρ # ≈ ρ̄ - @show α # ≈ ᾱ - @show β # ≈ β̄ - @show Q # ≈ Q̄ + @test ρ ≈ ρ̄ + @test α ≈ ᾱ + @test β ≈ β̄ + @test Q ≈ Q̄ - # :: Solutions :: + # Solutions x̂ = trunc.(Int, value.(x)) ŷ = objective_value(model) - @show x̂ # ∈ x̄ - @show ŷ # ≈ ȳ + @test x̂ ∈ x̄ + @test ŷ ≈ ȳ + + return nothing end end \ No newline at end of file diff --git a/test/examples/qba/qba2.jl b/test/examples/qba/qba2.jl index 67fea4f0..abd91b7c 100644 --- a/test/examples/qba/qba2.jl +++ b/test/examples/qba/qba2.jl @@ -1,5 +1,6 @@ function test_qba2() @testset "Illustrative Example" begin + # Problem Data Q̄ = [ -5 4 8 0 0 -3 2 0 @@ -12,6 +13,7 @@ function test_qba2() x̄ = Set{Vector{Int}}([[1, 0, 0, 1]]) ȳ = -11 + # Model model = Model(() -> ToQUBO.Optimizer(ExactSampler.Optimizer)) @variable(model, x[1:4], Bin) @@ -27,18 +29,20 @@ function test_qba2() optimize!(model) - # :: Reformulation :: - Q, α, β = ToQUBO.qubo(unsafe_backend(model), Matrix) + # Reformulation + Q, α, β = ToQUBO.qubo(model, Matrix) @test α ≈ ᾱ @test β ≈ β̄ @test Q ≈ Q̄ - # :: Solution :: + # Solution x̂ = value.(x) ŷ = objective_value(model) @test x̂ ∈ x̄ @test ŷ ≈ ȳ + + return nothing end end \ No newline at end of file diff --git a/test/examples/qba/qba3_1.jl b/test/examples/qba/qba3_1.jl index dd44e0ed..bc9fd981 100644 --- a/test/examples/qba/qba3_1.jl +++ b/test/examples/qba/qba3_1.jl @@ -8,11 +8,11 @@ function test_qba3_1() such that the subset sums are as close to each other as possible. =# - # :: Data :: + # Problem Data S = Int[25, 7, 13, 31, 42, 17, 21, 10] m = 8 - # :: Results :: + # Results Q̄ = [ -3525 350 650 1550 2100 850 1050 500 0 -1113 182 434 588 238 294 140 @@ -34,7 +34,7 @@ function test_qba3_1() ]) ȳ = -6889 - # :: Model :: + # Model model = Model(() -> ToQUBO.Optimizer(ExactSampler.Optimizer)) @variable(model, x[1:m], Bin) @@ -42,17 +42,19 @@ function test_qba3_1() optimize!(model) - Q, _, c = ToQUBO.qubo(unsafe_backend(model), Matrix) + Q, _, c = ToQUBO.qubo(model, Matrix) - # :: Reformulation :: + # Reformulation @test c ≈ c̄ @test Q ≈ 4Q̄ - # :: Solution :: + # Solution x̂ = trunc.(Int, value.(x)) ŷ = objective_value(model) @test x̂ ∈ x̄ @test ŷ ≈ 4ȳ + c̄ + + return nothing end end \ No newline at end of file diff --git a/test/examples/qba/qba3_2.jl b/test/examples/qba/qba3_2.jl index 66cad564..cd184ed3 100644 --- a/test/examples/qba/qba3_2.jl +++ b/test/examples/qba/qba3_2.jl @@ -18,7 +18,7 @@ function test_qba3_2() ⊻(x::VariableRef, y::VariableRef) = x + y - 2 * x * y - # :: Data :: + # Problem Data G = Dict{Tuple{Int,Int},Float64}( (1, 2) => 1.0, (1, 3) => 1.0, @@ -29,7 +29,7 @@ function test_qba3_2() ) m = 5 - # :: Results :: + # Results Q̄ = [ 2 -2 -2 0 0 0 2 0 -2 0 @@ -47,7 +47,7 @@ function test_qba3_2() ]) ȳ = 5 - # :: Model :: + # Model model = Model(() -> ToQUBO.Optimizer(ExactSampler.Optimizer)) @variable(model, x[1:m], Bin) @@ -55,17 +55,19 @@ function test_qba3_2() optimize!(model) - Q, _, c = ToQUBO.qubo(unsafe_backend(model), Matrix) + Q, _, c = ToQUBO.qubo(model, Matrix) - # :: Reformulation :: + # Reformulation @test c ≈ c̄ @test Q ≈ Q̄ - # :: Solutions :: + # Solutions x̂ = trunc.(Int, value.(x)) ŷ = objective_value(model) @test x̂ ∈ x̄ @test ŷ ≈ ȳ + + return nothing end end \ No newline at end of file diff --git a/test/examples/quadratic/quadratic.jl b/test/examples/quadratic/quadratic.jl index 8755532d..f2a15b12 100644 --- a/test/examples/quadratic/quadratic.jl +++ b/test/examples/quadratic/quadratic.jl @@ -1,7 +1,7 @@ -include("quadratic1.jl") +include("quadratic_1.jl") function test_quadratic() @testset "Quadratic Programs" verbose = true begin - test_quadratic1() + test_quadratic_1() end end \ No newline at end of file diff --git a/test/examples/quadratic/quadratic1.jl b/test/examples/quadratic/quadratic1.jl deleted file mode 100644 index 9770c955..00000000 --- a/test/examples/quadratic/quadratic1.jl +++ /dev/null @@ -1,75 +0,0 @@ -function test_quadratic1() - @testset "3 variables, 1 constraint" begin - # ~*~ Problem Data ~*~ # - n = 3 - A = [ - -1 2 2 - 2 -1 2 - 2 2 -1 - ] - b = 6 - - # Penalty Choice - ρ̄ = -16 - - # ~*~ Solution Data ~*~ # - Q̄ = [ - -209 740 740 32 64 128 64 -512 0 0 0 -1152 -512 0 -128 -256 -256 -256 -256 -128 - 0 -209 -412 -96 -192 -384 -192 512 -256 -512 -256 1152 0 -128 0 0 0 256 256 128 - 0 0 -209 -224 -448 -896 -448 0 256 512 256 1152 512 128 128 256 256 0 0 0 - 0 0 0 176 -64 -128 -64 0 0 0 0 0 0 128 128 0 0 0 0 128 - 0 0 0 0 320 -256 -128 0 256 0 0 0 0 0 0 256 0 256 0 0 - 0 0 0 0 0 512 -256 512 0 512 0 0 512 0 0 0 0 0 0 0 - 0 0 0 0 0 0 320 0 0 0 256 0 0 0 0 0 256 0 256 0 - 0 0 0 0 0 0 0 -512 0 0 0 0 0 0 0 0 0 0 0 0 - 0 0 0 0 0 0 0 0 -256 0 0 0 0 0 0 0 0 0 0 0 - 0 0 0 0 0 0 0 0 0 -512 0 0 0 0 0 0 0 0 0 0 - 0 0 0 0 0 0 0 0 0 0 -256 0 0 0 0 0 0 0 0 0 - 0 0 0 0 0 0 0 0 0 0 0 -1152 0 0 0 0 0 0 0 0 - 0 0 0 0 0 0 0 0 0 0 0 0 -512 0 0 0 0 0 0 0 - 0 0 0 0 0 0 0 0 0 0 0 0 0 -128 0 0 0 0 0 0 - 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -128 0 0 0 0 0 - 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -256 0 0 0 0 - 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -256 0 0 0 - 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -256 0 0 - 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -256 0 - 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -128 - ] - - ᾱ = 1 - β̄ = -576 - - x̄ = Set{Vector{Int}}([ - [0, 1, 1], - [1, 0, 1], - [1, 1, 0], - ]) - ȳ = 2 - - # ~*~ Model ~*~ # - model = Model(() -> ToQUBO.Optimizer(ExactSampler.Optimizer)) - - @variable(model, x[1:n], Bin) - @objective(model, Max, x'A * x) - @constraint(model, c1, x'A * x <= b) - - optimize!(model) - - # :: Reformulation :: - m = unsafe_backend(model) - ρ = MOI.get(m, ToQUBO.Penalty(), c1.index) - Q, α, β = ToQUBO.qubo(m, Matrix) - - @test ρ ≈ ρ̄ - @test α ≈ ᾱ - @test β ≈ β̄ - @test Q ≈ Q̄ broken = true - - # :: Solutions :: - x̂ = trunc.(Int, value.(x)) - ŷ = objective_value(model) - - @test x̂ ∈ x̄ - @test ŷ ≈ ȳ - end -end \ No newline at end of file diff --git a/test/examples/quadratic/quadratic_1.jl b/test/examples/quadratic/quadratic_1.jl new file mode 100644 index 00000000..5e45c6cc --- /dev/null +++ b/test/examples/quadratic/quadratic_1.jl @@ -0,0 +1,74 @@ +function test_quadratic_1() + @testset "3 variables, 1 constraint" begin + # Problem Data + n = 3 + A = [ + -1 2 2 + 2 -1 2 + 2 2 -1 + ] + b = 6 + + # Penalty Choice + ρ̄ = -16 + + # Solution + Q̄ = [ + -209 740 740 32 64 128 64 -1152 -128 -256 -512 -256 -128 -256 -512 -256 0 0 0 0 + 0 -209 -412 -96 -192 -384 -192 1152 128 256 512 256 0 0 0 0 -128 -256 -512 -256 + 0 0 -209 -224 -448 -896 -448 1152 0 0 0 0 128 256 512 256 128 256 512 256 + 0 0 0 176 -64 -128 -64 0 128 0 0 0 128 0 0 0 128 0 0 0 + 0 0 0 0 320 -256 -128 0 0 256 0 0 0 256 0 0 0 256 0 0 + 0 0 0 0 0 512 -256 0 0 0 512 0 0 0 512 0 0 0 512 0 + 0 0 0 0 0 0 320 0 0 0 0 256 0 0 0 256 0 0 0 256 + 0 0 0 0 0 0 0 -1152 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 -128 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 -256 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 -512 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 -256 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 -128 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 -256 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -512 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -256 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -128 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -256 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -512 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -256 + ] + + ᾱ = 1 + β̄ = -576 + + x̄ = Set{Vector{Int}}([[0, 1, 1], [1, 0, 1], [1, 1, 0]]) + ȳ = 2 + + # Model + model = Model(() -> ToQUBO.Optimizer(ExactSampler.Optimizer)) + + @variable(model, x[1:n], Bin) + @objective(model, Max, x'A * x) + @constraint(model, c1, x'A * x <= b) + + set_optimizer_attribute(model, TQA.StableQuadratization(), true) + + optimize!(model) + + # Reformulation + ρ = MOI.get(model, TQA.ConstraintEncodingPenalty(), c1) + Q, α, β = ToQUBO.qubo(model, Matrix) + + @test ρ ≈ ρ̄ + @test α ≈ ᾱ + @test β ≈ β̄ + @test Q ≈ Q̄ + + # Solutions + x̂ = trunc.(Int, value.(x)) + ŷ = objective_value(model) + + @test x̂ ∈ x̄ + @test ŷ ≈ ȳ + + return nothing + end +end \ No newline at end of file diff --git a/test/integration/attributes.jl b/test/integration/attributes.jl new file mode 100644 index 00000000..8f20b36a --- /dev/null +++ b/test/integration/attributes.jl @@ -0,0 +1,26 @@ +function test_attributes() + @testset "Attributes" begin + model = JuMP.Model(() -> ToQUBO.Optimizer(RandomSampler.Optimizer)) + + # MOI Attributes + @test MOI.get(model, MOI.NumberOfVariables()) == 0 + + @test MOI.get(model, MOI.TimeLimitSec()) |> isnothing + MOI.set(model, MOI.TimeLimitSec(), 1.0) + @test MOI.get(model, MOI.TimeLimitSec()) == 1.0 + + # ToQUBO Attributes + @test MOI.get(model, TQA.Quadratize()) === false + MOI.set(model, TQA.Quadratize(), true) + @test MOI.get(model, TQA.Quadratize()) === true + + # Solver Attributes + @test MOI.get(model, RandomSampler.RandomSeed()) |> isnothing + MOI.set(model, RandomSampler.RandomSeed(), 13) + @test MOI.get(model, RandomSampler.RandomSeed()) == 13 + + @test MOI.get(model, RandomSampler.NumberOfReads()) == 1_000 + MOI.set(model, RandomSampler.NumberOfReads(), 13) + @test MOI.get(model, RandomSampler.NumberOfReads()) == 13 + end +end \ No newline at end of file diff --git a/test/integration/integration.jl b/test/integration/integration.jl index 2dcf9db6..85856c62 100644 --- a/test/integration/integration.jl +++ b/test/integration/integration.jl @@ -1,3 +1,7 @@ -function test_integration() +include("attributes.jl") +function test_integration() + @testset "Integration" verbose = true begin + test_attributes() + end end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 49e34bc9..cba91989 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,35 +1,24 @@ using Test -# -*- Imports: JuMP + MOI -*- using JuMP const MOIU = MOI.Utilities const VI = MOI.VariableIndex const CI = MOI.ConstraintIndex -# -*- Imports: JuMP + MOI -*- -import MutableArithmetics -const MA = MutableArithmetics - -# -*- Imports -*- using ToQUBO: ToQUBO, PBO using Anneal using LinearAlgebra using TOML -# -*- Test Assets -*- # -include("assets/assets.jl") +const TQA = ToQUBO.Attributes -# -*- Unit Tests -*- # +include("assets/assets.jl") include("unit/unit.jl") - -# -*- Integration Tests -*- # include("integration/integration.jl") - -# -*- Examples -*- # include("examples/examples.jl") function main() - @testset ":: -*- :: ~*~ :: ToQUBO.jl :: ~*~ :: -*- ::" verbose = true begin + @testset ":: :: :: ToQUBO.jl :: :: ::" verbose = true begin test_unit() test_integration() test_examples() diff --git a/test/unit/compiler/constraints.jl b/test/unit/compiler/constraints.jl index 3444a6c0..512917ce 100644 --- a/test/unit/compiler/constraints.jl +++ b/test/unit/compiler/constraints.jl @@ -7,7 +7,7 @@ function test_compiler_constraints_quadratic() ] b = 6.0 - model = ToQUBO.VirtualQUBOModel{Float64}() + model = ToQUBO.VirtualModel{Float64}() arch = ToQUBO.GenericArchitecture() x, _ = MOI.add_constrained_variables(model.source_model, repeat([MOI.ZeroOne()], n)) @@ -21,6 +21,20 @@ function test_compiler_constraints_quadratic() s = MOI.EqualTo{Float64}(b) g = ToQUBO.toqubo_constraint(model, f, s, arch) + h = ToQUBO.PBO.PBF{VI,Float64}() + + ToQUBO.toqubo_parse!(model, h, f, s, arch) + + @test h == PBO.PBF{VI,Float64}( + -6.0, + x[1] => -0.5, + x[2] => -0.5, + x[3] => -0.5, + [x[1], x[2]] => 4.0, + [x[1], x[3]] => 4.0, + [x[2], x[3]] => 4.0, + ) + @test g == PBO.PBF{VI,Float64}( 144.0, x[1] => 25.0, diff --git a/test/unit/lib/pbo.jl b/test/unit/lib/pbo.jl index 0447013f..57922b20 100644 --- a/test/unit/lib/pbo.jl +++ b/test/unit/lib/pbo.jl @@ -1,7 +1,4 @@ function test_pbo() - S = Symbol - T = Float64 - @testset "Pseudo-Boolean Optimization" verbose = true begin @testset "Constructors" begin for (x, y) in Assets.PBF_CONSTRUCTOR_LIST @@ -45,7 +42,7 @@ end # function test_pbo() # @testset "PBO" verbose = true begin -# # -*- Definitions -*- +# # Definitions # S = Symbol # T = Float64 diff --git a/test/unit/lib/pbo_ma.jl b/test/unit/lib/pbo_ma.jl deleted file mode 100644 index cd684c6e..00000000 --- a/test/unit/lib/pbo_ma.jl +++ /dev/null @@ -1,66 +0,0 @@ -function test_pbo_ma() - @testset "PBO MutableArithmetics" verbose = true begin - # -*- Definitions -*- - S = Symbol - T = Float64 - - # :: Canonical Constructor :: - f = PBO.PBF{S,T}( - Dict{Set{S},T}( - Set{S}() => 0.5, - Set{S}([:x]) => 1.0, - Set{S}([:y]) => 1.0, - Set{S}([:z]) => 1.0, - Set{S}([:x, :y]) => -2.0, - Set{S}([:x, :z]) => -2.0, - Set{S}([:y, :z]) => -2.0, - Set{S}([:x, :y, :z]) => 3.0, - ), - ) - g = PBO.PBF{S,T}(Dict{Set{S},T}(Set{S}() => 1.0)) - h = PBO.PBF{S,T}( - Dict{Set{S},T}( - Set{S}([:x]) => 1.0, - Set{S}([:y]) => 1.0, - Set{S}([:z]) => 1.0, - Set{S}([:w]) => 1.0, - Set{S}() => 1.0, - ), - ) - p = PBO.PBF{S,T}( - Dict{Set{S},T}(Set{S}() => 0.5, Set{S}([:x]) => 1.0, Set{S}([:x, :y]) => -2.0), - ) - q = PBO.PBF{S,T}( - Dict{Set{S},T}(Set{S}() => 0.5, Set{S}([:y]) => 1.0, Set{S}([:x, :y]) => 2.0), - ) - r = PBO.PBF{S,T}(Set{S}() => 1.0, Set{S}([:z]) => -1.0) - s = PBO.PBF{S,T}(Set{S}() => 0.0, Set{S}([:x, :y, :z]) => 3.0) - - @testset "Arithmetic" verbose = true begin - MA.Test.@test_rewrite(p + q) - MA.Test.@test_rewrite(p + q + r) - MA.Test.@test_rewrite(s + 3.0) - MA.Test.@test_rewrite(p - q) - MA.Test.@test_rewrite(s - 3.0) - MA.Test.@test_rewrite(0.0 - f) - MA.Test.@test_rewrite(3.0 - s) - MA.Test.@test_rewrite(p + q - r) - MA.Test.@test_rewrite(p * q) - MA.Test.@test_rewrite(p + 2.0*q) - MA.Test.@test_rewrite(p - 2.0*q) - MA.Test.@test_rewrite(p*2.0) - MA.Test.@test_rewrite(f*0.0) - MA.Test.@test_rewrite(f*-2.0) - MA.Test.@test_rewrite(p / 2.0) - MA.Test.@test_rewrite(p^0) - MA.Test.@test_rewrite(q^2) - MA.Test.@test_rewrite(r^2) - MA.Test.@test_rewrite(s^2) - MA.Test.@test_rewrite(r^3) - MA.Test.@test_rewrite(s^3) - MA.Test.@test_rewrite(r^4) - - end - - end -end \ No newline at end of file diff --git a/test/unit/lib/virtual.jl b/test/unit/lib/virtual.jl index a2cf9b0f..7962b58e 100644 --- a/test/unit/lib/virtual.jl +++ b/test/unit/lib/virtual.jl @@ -1,48 +1,20 @@ -struct VirtualModel{T} <: ToQUBO.AbstractVirtualModel{T} - source_model::MOIU.Model{T} - target_model::MOIU.Model{T} - variables::Vector{ToQUBO.VirtualVariable{T}} - source::Dict{VI,ToQUBO.VirtualVariable{T}} - target::Dict{VI,ToQUBO.VirtualVariable{T}} - - function VirtualModel{T}(; kws...) where {T} - new{T}( - MOIU.Model{T}(), - MOIU.Model{T}(), - ToQUBO.VirtualVariable{T}[], - Dict{VI,ToQUBO.VirtualVariable{T}}(), - Dict{VI,ToQUBO.VirtualVariable{T}}(), - ) - end - - VirtualModel(; kws...) = VirtualModel{Float64}(; kws...) -end - -MOI.get(model::VirtualModel, ::ToQUBO.SourceModel) = model.source_model -MOI.get(model::VirtualModel, ::ToQUBO.TargetModel) = model.target_model -MOI.get(model::VirtualModel, ::ToQUBO.Variables) = model.variables -MOI.get(model::VirtualModel, ::ToQUBO.Source, x::VI) = model.source[x] -MOI.set(model::VirtualModel, ::ToQUBO.Source, x::VI, v::ToQUBO.VV) = (model.source[x] = v) -MOI.get(model::VirtualModel, ::ToQUBO.Target, y::VI) = model.target[y] -MOI.set(model::VirtualModel, ::ToQUBO.Target, y::VI, v::ToQUBO.VV) = (model.target[y] = v) - function test_virtual() @testset "VirtualMapping" verbose = true begin @testset "Virtual Model" begin - model = VirtualModel() + model = ToQUBO.VirtualModel() @test MOI.is_empty(model) end @testset "Encodings" verbose = true begin @testset "Linear" begin - model = VirtualModel() + model = ToQUBO.VirtualModel() - x = MOI.add_variable(MOI.get(model, ToQUBO.SourceModel())) + x = MOI.add_variable(model.source_model) γ = [1.0, 2.0, 3.0] α = 1.0 - v = ToQUBO.encode!(ToQUBO.Linear(), model, x, γ, α) + v = ToQUBO.encode!(model, ToQUBO.Linear(), x, γ, α) y = ToQUBO.target(v) @test length(y) == 3 @@ -56,17 +28,17 @@ function test_virtual() ) @test isnothing(ToQUBO.penaltyfn(v)) - @test MOI.get(model, ToQUBO.Variables()) == [v] - @test MOI.get(model, ToQUBO.Source(), ToQUBO.source(v)) == v - @test MOI.get.(model, ToQUBO.Target(), ToQUBO.target(v)) == [v, v, v] + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v] end @testset "Mirror" begin - model = VirtualModel() + model = ToQUBO.VirtualModel() - x = MOI.add_variable(MOI.get(model, ToQUBO.SourceModel())) + x = MOI.add_variable(model.source_model) - v = ToQUBO.encode!(ToQUBO.Mirror(), model, x) + v = ToQUBO.encode!(model, ToQUBO.Mirror(), x) y = ToQUBO.target(v) @test length(y) == 1 @@ -75,18 +47,18 @@ function test_virtual() @test ToQUBO.expansion(v) == PBO.PBF{VI,Float64}(y[1]) @test isnothing(ToQUBO.penaltyfn(v)) - @test MOI.get(model, ToQUBO.Variables()) == [v] - @test MOI.get(model, ToQUBO.Source(), ToQUBO.source(v)) == v - @test MOI.get.(model, ToQUBO.Target(), ToQUBO.target(v)) == [v] + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v] end @testset "Unary ℤ" begin - model = VirtualModel() + model = ToQUBO.VirtualModel() - x = MOI.add_variable(MOI.get(model, ToQUBO.SourceModel())) + x = MOI.add_variable(model.source_model) a, b = (-2.0, 2.0) - v = ToQUBO.encode!(ToQUBO.Unary(), model, x, a, b) + v = ToQUBO.encode!(model, ToQUBO.Unary(), x, a, b) y = ToQUBO.target(v) @test length(y) == 4 @@ -101,19 +73,19 @@ function test_virtual() ) @test isnothing(ToQUBO.penaltyfn(v)) - @test MOI.get(model, ToQUBO.Variables()) == [v] - @test MOI.get(model, ToQUBO.Source(), ToQUBO.source(v)) == v - @test MOI.get.(model, ToQUBO.Target(), ToQUBO.target(v)) == [v, v, v, v] + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v, v] end @testset "Unary ℝ" begin - model = VirtualModel() + model = ToQUBO.VirtualModel() - x = MOI.add_variable(MOI.get(model, ToQUBO.SourceModel())) + x = MOI.add_variable(model.source_model) a, b = (-2.0, 2.0) n = 4 - v = ToQUBO.encode!(ToQUBO.Unary(), model, x, a, b, n) + v = ToQUBO.encode!(model, ToQUBO.Unary(), x, a, b, n) y = ToQUBO.target(v) @test length(y) == n @@ -128,18 +100,18 @@ function test_virtual() ) @test isnothing(ToQUBO.penaltyfn(v)) - @test MOI.get(model, ToQUBO.Variables()) == [v] - @test MOI.get(model, ToQUBO.Source(), ToQUBO.source(v)) == v - @test MOI.get.(model, ToQUBO.Target(), ToQUBO.target(v)) == [v, v, v, v] + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v, v] end @testset "Binary ℤ" begin - model = VirtualModel() + model = ToQUBO.VirtualModel() - x = MOI.add_variable(MOI.get(model, ToQUBO.SourceModel())) + x = MOI.add_variable(model.source_model) a, b = (-2.0, 2.0) - v = ToQUBO.encode!(ToQUBO.Binary(), model, x, a, b) + v = ToQUBO.encode!(model, ToQUBO.Binary(), x, a, b) y = ToQUBO.target(v) @test length(y) == 3 @@ -153,19 +125,19 @@ function test_virtual() ) @test isnothing(ToQUBO.penaltyfn(v)) - @test MOI.get(model, ToQUBO.Variables()) == [v] - @test MOI.get(model, ToQUBO.Source(), ToQUBO.source(v)) == v - @test MOI.get.(model, ToQUBO.Target(), ToQUBO.target(v)) == [v, v, v] + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v] end @testset "Binary ℝ" begin - model = VirtualModel() + model = ToQUBO.VirtualModel() - x = MOI.add_variable(MOI.get(model, ToQUBO.SourceModel())) + x = MOI.add_variable(model.source_model) a, b = (-2.0, 2.0) n = 3 - v = ToQUBO.encode!(ToQUBO.Binary(), model, x, a, b, n) + v = ToQUBO.encode!(model, ToQUBO.Binary(), x, a, b, n) y = ToQUBO.target(v) @test length(y) == n @@ -179,18 +151,18 @@ function test_virtual() ) @test isnothing(ToQUBO.penaltyfn(v)) - @test MOI.get(model, ToQUBO.Variables()) == [v] - @test MOI.get(model, ToQUBO.Source(), ToQUBO.source(v)) == v - @test MOI.get.(model, ToQUBO.Target(), ToQUBO.target(v)) == [v, v, v] + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v] end @testset "Arithmetic ℤ" begin - model = VirtualModel() + model = ToQUBO.VirtualModel() - x = MOI.add_variable(MOI.get(model, ToQUBO.SourceModel())) + x = MOI.add_variable(model.source_model) a, b = (-2.0, 2.0) - v = ToQUBO.encode!(ToQUBO.Arithmetic(), model, x, a, b) + v = ToQUBO.encode!(model, ToQUBO.Arithmetic(), x, a, b) y = ToQUBO.target(v) @test length(y) == 3 @@ -204,19 +176,19 @@ function test_virtual() ) @test isnothing(ToQUBO.penaltyfn(v)) - @test MOI.get(model, ToQUBO.Variables()) == [v] - @test MOI.get(model, ToQUBO.Source(), ToQUBO.source(v)) == v - @test MOI.get.(model, ToQUBO.Target(), ToQUBO.target(v)) == [v, v, v] + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v] end @testset "Arithmetic ℝ" begin - model = VirtualModel() + model = ToQUBO.VirtualModel() - x = MOI.add_variable(MOI.get(model, ToQUBO.SourceModel())) + x = MOI.add_variable(model.source_model) a, b = (-2.0, 2.0) n = 3 - v = ToQUBO.encode!(ToQUBO.Arithmetic(), model, x, a, b, n) + v = ToQUBO.encode!(model, ToQUBO.Arithmetic(), x, a, b, n) y = ToQUBO.target(v) @test length(y) == n @@ -230,18 +202,18 @@ function test_virtual() ) @test isnothing(ToQUBO.penaltyfn(v)) - @test MOI.get(model, ToQUBO.Variables()) == [v] - @test MOI.get(model, ToQUBO.Source(), ToQUBO.source(v)) == v - @test MOI.get.(model, ToQUBO.Target(), ToQUBO.target(v)) == [v, v, v] + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v] end @testset "One Hot" begin - model = VirtualModel() + model = ToQUBO.VirtualModel() - x = MOI.add_variable(MOI.get(model, ToQUBO.SourceModel())) + x = MOI.add_variable(model.source_model) γ = [-1.0, -0.5, 0.0, 0.5, 1.0] - v = ToQUBO.encode!(ToQUBO.OneHot(), model, x, γ) + v = ToQUBO.encode!(model, ToQUBO.OneHot(), x, γ) y = ToQUBO.target(v) @test length(y) == length(γ) @@ -255,18 +227,18 @@ function test_virtual() ) @test ToQUBO.penaltyfn(v) ≈ (PBO.PBF{VI,Float64}(-1.0, y...)^2) - @test MOI.get(model, ToQUBO.Variables()) == [v] - @test MOI.get(model, ToQUBO.Source(), ToQUBO.source(v)) == v - @test MOI.get.(model, ToQUBO.Target(), ToQUBO.target(v)) == [v, v, v, v, v] + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v, v, v] end @testset "One Hot ℤ" begin - model = VirtualModel() + model = ToQUBO.VirtualModel() - x = MOI.add_variable(MOI.get(model, ToQUBO.SourceModel())) + x = MOI.add_variable(model.source_model) a, b = (-2.0, 2.0) - v = ToQUBO.encode!(ToQUBO.OneHot(), model, x, a, b) + v = ToQUBO.encode!(model, ToQUBO.OneHot(), x, a, b) y = ToQUBO.target(v) @test length(y) == 5 @@ -280,18 +252,18 @@ function test_virtual() ) @test ToQUBO.penaltyfn(v) ≈ (PBO.PBF{VI,Float64}(-1.0, y...)^2) - @test MOI.get(model, ToQUBO.Variables()) == [v] - @test MOI.get(model, ToQUBO.Source(), ToQUBO.source(v)) == v - @test MOI.get.(model, ToQUBO.Target(), ToQUBO.target(v)) == [v, v, v, v, v] + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v, v, v] end @testset "One Hot ℝ" begin - model = VirtualModel() + model = ToQUBO.VirtualModel() - x = MOI.add_variable(MOI.get(model, ToQUBO.SourceModel())) + x = MOI.add_variable(model.source_model) a, b = (-2.0, 2.0) n = 5 - v = ToQUBO.encode!(ToQUBO.OneHot(), model, x, a, b, n) + v = ToQUBO.encode!(model, ToQUBO.OneHot(), x, a, b, n) y = ToQUBO.target(v) @test length(y) == n @@ -305,17 +277,17 @@ function test_virtual() ) @test ToQUBO.penaltyfn(v) ≈ (PBO.PBF{VI,Float64}(-1.0, y...)^2) - @test MOI.get(model, ToQUBO.Variables()) == [v] - @test MOI.get(model, ToQUBO.Source(), ToQUBO.source(v)) == v - @test MOI.get.(model, ToQUBO.Target(), ToQUBO.target(v)) == [v, v, v, v, v] + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v, v, v] end @testset "Domain Wall ℤ" begin - model = VirtualModel() + model = ToQUBO.VirtualModel() - x = MOI.add_variable(MOI.get(model, ToQUBO.SourceModel())) + x = MOI.add_variable(model.source_model) a, b = (-2.0, 2.0) - v = ToQUBO.encode!(ToQUBO.DomainWall(), model, x, a, b) + v = ToQUBO.encode!(model, ToQUBO.DomainWall(), x, a, b) y = ToQUBO.target(v) @test length(y) == 4 @@ -336,19 +308,19 @@ function test_virtual() [y[3], y[4]] => -2.0, ) - @test MOI.get(model, ToQUBO.Variables()) == [v] - @test MOI.get(model, ToQUBO.Source(), ToQUBO.source(v)) == v - @test MOI.get.(model, ToQUBO.Target(), ToQUBO.target(v)) == [v, v, v, v] + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v, v] end @testset "Domain Wall ℝ" begin - model = VirtualModel() + model = ToQUBO.VirtualModel() - x = MOI.add_variable(MOI.get(model, ToQUBO.SourceModel())) + x = MOI.add_variable(model.source_model) a, b = (-2.0, 2.0) n = 5 - v = ToQUBO.encode!(ToQUBO.DomainWall(), model, x, a, b, n) + v = ToQUBO.encode!(model, ToQUBO.DomainWall(), x, a, b, n) y = ToQUBO.target(v) @test length(y) == n - 1 @@ -369,9 +341,92 @@ function test_virtual() [y[3], y[4]] => -2.0, ) - @test MOI.get(model, ToQUBO.Variables()) == [v] - @test MOI.get(model, ToQUBO.Source(), ToQUBO.source(v)) == v - @test MOI.get.(model, ToQUBO.Target(), ToQUBO.target(v)) == [v, v, v, v] + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v, v] + end + + @testset "Bounded(Unary) ℤ" begin + model = ToQUBO.VirtualModel() + + x = MOI.add_variable(model.source_model) + a, b = (-10.0, 10.0) + + v = ToQUBO.encode!(model, ToQUBO.Bounded{ToQUBO.Unary}(5.0), x, a, b) + y = ToQUBO.target(v) + + @test length(y) == 8 + @test ToQUBO.source(v) == x + @test ToQUBO.expansion(v) == PBO.PBF{VI,Float64}( + y[1] => 1.0, + y[2] => 1.0, + y[3] => 1.0, + y[4] => 1.0, + y[5] => 5.0, + y[6] => 5.0, + y[7] => 5.0, + y[8] => 1.0, + nothing => a, + ) + @test isnothing(ToQUBO.penaltyfn(v)) + + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v, v, v, v, v, v] + end + + @testset "Bounded(Binary) ℤ" begin + model = ToQUBO.VirtualModel() + + x = MOI.add_variable(model.source_model) + a, b = (-10.0, 10.0) + + v = ToQUBO.encode!(model, ToQUBO.Bounded{ToQUBO.Binary}(5.0), x, a, b) + y = ToQUBO.target(v) + + @test length(y) == 6 + @test ToQUBO.source(v) == x + @test ToQUBO.expansion(v) == PBO.PBF{VI,Float64}( + y[1] => 1.0, + y[2] => 2.0, + y[3] => 4.0, + y[4] => 5.0, + y[5] => 5.0, + y[6] => 3.0, + nothing => a, + ) + @test isnothing(ToQUBO.penaltyfn(v)) + + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v, v, v, v] + end + + @testset "Bounded(Arithmetic) ℤ" begin + model = ToQUBO.VirtualModel() + + x = MOI.add_variable(model.source_model) + a, b = (-10.0, 10.0) + + v = ToQUBO.encode!(model, ToQUBO.Bounded{ToQUBO.Arithmetic}(5.0), x, a, b) + y = ToQUBO.target(v) + + @test length(y) == 6 + @test ToQUBO.source(v) == x + @test ToQUBO.expansion(v) == PBO.PBF{VI,Float64}( + y[1] => 1.0, + y[2] => 2.0, + y[3] => 3.0, + y[4] => 4.0, + y[5] => 5.0, + y[6] => 5.0, + nothing => a, + ) + @test isnothing(ToQUBO.penaltyfn(v)) + + @test model.variables == [v] + @test model.source[ToQUBO.source(v)] == (v) + @test [model.target[y] for y in ToQUBO.target(v)] == [v, v, v, v, v, v] end end end