diff --git a/docs/src/examples.md b/docs/src/examples.md index e3ec4c3..6fa8107 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -141,6 +141,38 @@ selected = assets[findall(>(0.5), value.(x))] (true, "Wind, Battery") ``` +### Selecting the PEPS Backend from JuMP + +The optional SpinGlassPEPS backend is selected explicitly and requires +structured topology metadata. It is experimental and intended for structured +quasi-two-dimensional QUBOs; arbitrary dense QUBOs should remain on the default +DMRG backend. + +```julia +using JuMP, TenSolver + +m, n = 2, 2 +model = Model(TenSolver.Optimizer) +set_attribute(model, "backend", :peps) +set_attribute(model, "peps_layout", :square) +set_attribute(model, "peps_topology", (m, n)) +set_attribute(model, "peps_beta", 2.0) +set_attribute(model, "peps_bond_dim", 8) +set_attribute(model, "peps_max_states", 256) +set_attribute(model, "peps_cutoff_prob", 0.0) +set_attribute(model, "peps_strategy", :svd) +set_attribute(model, "peps_transformations", :identity) + +@variable(model, x[1:(m * n)], Bin) +# Add a structured objective whose variable order matches the square grid. +@objective(model, Min, -sum(x)) +optimize!(model) +``` + +If the optional SpinGlass component packages are not available, this backend +errors clearly. If `"backend"` is left unset, or set to `:dmrg`, TenSolver uses +the existing DMRG path. + ## Controlling Solver Parameters You can control various solver parameters for better performance: diff --git a/docs/src/spinglasspeps_integration.md b/docs/src/spinglasspeps_integration.md index a8169d1..cf4792e 100644 --- a/docs/src/spinglasspeps_integration.md +++ b/docs/src/spinglasspeps_integration.md @@ -207,8 +207,30 @@ Ising graph, clusters it with `super_square_lattice`, constructs the Potts Hamiltonian, runs `MpsContractor` plus `low_energy_spectrum`, and decodes retained states back to TenSolver Boolean vectors. -Later PRs should add QUBODrivers/JuMP raw optimizer attributes for backend and -PEPS parameters. +The QUBODrivers/JuMP optimizer interface can also select the backend through +raw optimizer attributes. DMRG remains the default: + +```julia +set_attribute(model, "backend", :dmrg) +``` + +The PEPS path is selected explicitly and requires topology metadata: + +```julia +set_attribute(model, "backend", :peps) +set_attribute(model, "peps_layout", :square) +set_attribute(model, "peps_topology", (m, n)) +set_attribute(model, "peps_beta", 2.0) +set_attribute(model, "peps_bond_dim", 8) +set_attribute(model, "peps_max_states", 256) +set_attribute(model, "peps_cutoff_prob", 0.0) +set_attribute(model, "peps_strategy", :svd) +``` + +PEPS runs add namespaced data under the returned QUBOTools `SampleSet` +metadata, including the selected backend, topology/layout, contraction/search +parameters, candidate-state count, effective time, selected transformation, and +largest discarded probability when available. Any PEPS selection API must validate that the problem includes enough topology metadata for the structured backend. If the topology is missing or unsupported, diff --git a/src/TenSolver.jl b/src/TenSolver.jl index 28f6da8..5cedd9e 100644 --- a/src/TenSolver.jl +++ b/src/TenSolver.jl @@ -47,12 +47,167 @@ QUBODrivers.@setup Optimizer begin EigsolveTol["eigsolve_tol"] :: Float64 = 1e-14 Preprocess["preprocess"] :: Bool = false Verbosity["verbosity"] :: Int = 1 + # Backend selection + "backend" :: Union{Symbol, String} = :dmrg + # PEPS backend keywords + "peps_topology" :: Any = nothing + "peps_layout" :: Union{Symbol, String} = :square + "peps_beta" :: Float64 = 2.0 + "peps_bond_dim" :: Int = 8 + "peps_max_states" :: Int = 256 + "peps_cutoff_prob" :: Float64 = 0.0 + "peps_onGPU" :: Bool = false + "peps_strategy" :: Union{Symbol, String} = :auto + "peps_num_sweeps" :: Int = 1 + "peps_transformations" :: Any = :all + "peps_local_dimension" :: Any = nothing end end QUBODrivers.honors_final_reads(::Type{<:Optimizer}) = true QUBODrivers.enforces_time_limit(::Type{<:Optimizer}) = true +function _optimizer_symbol(value::Symbol, attr::AbstractString) + return Symbol(lowercase(String(value))) +end + +function _optimizer_symbol(value::AbstractString, attr::AbstractString) + return Symbol(lowercase(strip(value))) +end + +function _optimizer_symbol(value, attr::AbstractString) + throw(ArgumentError("Optimizer attribute `$attr` must be a Symbol or String. Got $(repr(value)).")) +end + +function _optimizer_backend(get) + backend = _optimizer_symbol(get("backend"), "backend") + backend === :dmrg && return :dmrg + backend === :peps && return _optimizer_peps_backend(get) + throw(ArgumentError("Unsupported optimizer backend $(repr(get("backend"))). Use :dmrg or :peps.")) +end + +function _peps_topology_tuple(topology) + topology === nothing && + throw(ArgumentError("PEPS backend requires `peps_topology`, for example `(m, n)` for a square or king grid.")) + + if topology isa Tuple + return topology + elseif topology isa AbstractVector + return Tuple(topology) + else + throw(ArgumentError("`peps_topology` must be a tuple/vector such as `(m, n)` or `(m, n, spins_per_site)`. Got $(repr(topology)).")) + end +end + +function _peps_topology(layout, topology) + topology isa AbstractStructuredTopology && return topology + + dims = _peps_topology_tuple(topology) + if !(length(dims) in (2, 3)) + throw(ArgumentError("`peps_topology` must have 2 or 3 entries. Got $(repr(topology)).")) + end + + layout = _optimizer_symbol(layout, "peps_layout") + layout === :square && return SquareGrid(dims...) + layout === :king && return KingGrid(dims...) + throw(ArgumentError("Unsupported `peps_layout` $(repr(layout)). Use :square or :king.")) +end + +function _peps_local_dimension(local_dimension) + local_dimension === nothing && return nothing + local_dimension isa Integer && return Int(local_dimension) + throw(ArgumentError("Only integer `peps_local_dimension` values are currently supported as local dimension limits. Got $(repr(local_dimension)).")) +end + +function _optimizer_peps_backend(get) + return PEPSBackend( + _peps_topology(get("peps_layout"), get("peps_topology")); + beta = get("peps_beta"), + bond_dim = get("peps_bond_dim"), + max_states = get("peps_max_states"), + cutoff_prob = get("peps_cutoff_prob"), + onGPU = get("peps_onGPU"), + contraction = _optimizer_symbol(get("peps_strategy"), "peps_strategy"), + num_sweeps = get("peps_num_sweeps"), + transformations = get("peps_transformations"), + local_dimension = _peps_local_dimension(get("peps_local_dimension")), + ) +end + +function _qubo_samples(::Type{T}, psi::Solution, l, Q, a, b, num_reads::Integer) where {T} + samples = Vector{QUBOTools.Sample{T,Int}}(undef, num_reads) + for i in 1:num_reads + x = sample(psi) + E = QUBOTools.value(x, l, Q, a, b) + + samples[i] = QUBOTools.Sample{T,Int}(x, E) + end + + return samples +end + +function _peps_read_counts(psi::PEPSSolution, num_reads::Integer) + num_reads >= 0 || throw(ArgumentError("num_reads must be nonnegative. Got $num_reads.")) + + states = psi.states + isempty(states) && throw(ArgumentError("Cannot build QUBOTools samples from an empty PEPS solution.")) + + counts = zeros(Int, length(states)) + num_reads == 0 && return counts + + probabilities = psi.probabilities + if isempty(probabilities) + counts[begin] = Int(num_reads) + return counts + end + + length(probabilities) == length(states) || + throw(ArgumentError("PEPS probabilities length must match states length. Got $(length(probabilities)) probabilities for $(length(states)) states.")) + any(p -> p < 0, probabilities) && + throw(ArgumentError("PEPS probabilities must be nonnegative. Got $(repr(probabilities)).")) + + total = sum(probabilities) + total > 0 || throw(ArgumentError("PEPS probabilities must have positive total weight. Got $(repr(probabilities)).")) + + weights = (Float64.(probabilities) ./ Float64(total)) .* Int(num_reads) + counts .= floor.(Int, weights) + remaining = Int(num_reads) - sum(counts) + + if remaining > 0 + order = sortperm(collect(eachindex(weights)); by = i -> (weights[i] - counts[i], -i), rev = true) + for i in Iterators.take(order, remaining) + counts[i] += 1 + end + end + + return counts +end + +function _qubo_samples(::Type{T}, psi::PEPSSolution, l, Q, a, b, num_reads::Integer) where {T} + counts = _peps_read_counts(psi, num_reads) + samples = QUBOTools.Sample{T,Int}[] + sizehint!(samples, count(>(0), counts)) + + for (x, reads) in zip(psi.states, counts) + reads == 0 && continue + E = QUBOTools.value(x, l, Q, a, b) + push!(samples, QUBOTools.Sample{T,Int}(copy(x), E, reads)) + end + + return samples +end + +function _add_backend_metadata!(metadata::Dict{String,Any}, psi::PEPSSolution) + peps = copy(psi.metadata) + peps["candidate_states"] = length(psi.states) + peps["effective_time"] = metadata["time"]["effective"] + tensolver = get!(metadata, "tensolver", Dict{String,Any}()) + tensolver["peps"] = peps + return metadata +end + +_add_backend_metadata!(metadata::Dict{String,Any}, psi) = metadata + function QUBODrivers.sample(sampler::Optimizer{T}) where {T} # ~ Manage Attributes ~ # get(attr) = MOI.get(sampler, MOI.RawOptimizerAttribute(attr)) @@ -72,30 +227,34 @@ function QUBODrivers.sample(sampler::Optimizer{T}) where {T} n, l, Q, a, b = QUBOTools.qubo(sampler, :sparse; sense = :min) # min_x a*(x'Qx + l'x + b) # s.t. x in {0, 1}^n - results = @timed minimize(Q, l, b; - cutoff = get("cutoff"), - vtol = get("vtol"), - iterations = get("iterations"), - time_limit, - maxdim = get("maxdim"), - mindim = get("mindim"), - noise = get("noise"), - device = get("device"), - verbosity, - eigsolve_krylovdim = get("eigsolve_krylovdim"), - eigsolve_tol = get("eigsolve_tol"), - eigsolve_maxiter = get("eigsolve_maxiter"), - ) + backend = _optimizer_backend(get) + results = if backend isa PEPSBackend + @timed minimize(Q, l, b; + backend, + cutoff = get("cutoff"), + verbosity, + ) + else + @timed minimize(Q, l, b; + backend, + cutoff = get("cutoff"), + vtol = get("vtol"), + iterations = get("iterations"), + time_limit, + maxdim = get("maxdim"), + mindim = get("mindim"), + noise = get("noise"), + device = get("device"), + verbosity, + eigsolve_krylovdim = get("eigsolve_krylovdim"), + eigsolve_tol = get("eigsolve_tol"), + eigsolve_maxiter = get("eigsolve_maxiter"), + ) + end _, psi = results.value # ~ Samples and Output ~ # - samples = Vector{QUBOTools.Sample{T,Int}}(undef, final_num_reads) - for i in 1:final_num_reads - x = sample(psi) - E = QUBOTools.value(x, l, Q, a, b) - - samples[i] = QUBOTools.Sample{T,Int}(x, E) - end + samples = _qubo_samples(T, psi, l, Q, a, b, final_num_reads) # ~ Metadata ~ # metadata = _tensolver_metadata( @@ -109,6 +268,7 @@ function QUBODrivers.sample(sampler::Optimizer{T}) where {T} vtol = get("vtol"), maxdim = get("maxdim"), ) + _add_backend_metadata!(metadata, psi) return QUBOTools.SampleSet{T}(samples, metadata; sense = :min, domain = :bool) end @@ -161,6 +321,45 @@ function _tensolver_metadata( return metadata end +function _tensolver_metadata( + solution::PEPSSolution; + effective_time::Real, + num_reads::Integer, + final_num_reads::Integer, + time_limit::Real, + iterations::Integer, + cutoff::Real, + vtol::Real, + maxdim, +) + algorithm_name = get(solution.metadata, "backend", "SpinGlassPEPS") + metadata = QUBODrivers._sampler_metadata( + origin = "TenSolver.jl", + algorithm_name = algorithm_name, + backend_name = "TenSolver", + backend_version = __VERSION__, + execution_mode = "tensor_network_peps", + optimizer_iterations = 1, + optimizer_evaluations = length(solution.states), + number_of_reads = num_reads, + final_number_of_reads = final_num_reads, + status = "locally_solved", + termination_status = MOI.LOCALLY_SOLVED, + ) + metadata["time"] = Dict{String,Any}("effective" => effective_time) + metadata["tensolver"] = Dict{String,Any}( + "parameters" => Dict{String,Any}( + "cutoff" => cutoff, + "vtol" => vtol, + "maxdim" => maxdim isa AbstractVector ? copy(maxdim) : maxdim, + "iterations" => iterations, + "time_limit" => time_limit, + ), + ) + + return metadata +end + function _tensolver_status(solution::Solution; iterations::Integer, time_limit::Real) elapsed_time = isempty(solution.elapsed_times) ? 0.0 : last(solution.elapsed_times) if length(solution.energies) >= iterations diff --git a/test/Project.toml b/test/Project.toml index ce8f394..14d6300 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -9,6 +9,7 @@ JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" QUBODrivers = "a3f166f7-2cd3-47b6-9e1e-6fbfe0449eb0" +QUBOTools = "60eb5b62-0a39-4ddc-84c5-97d2adff9319" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" diff --git a/test/jump.jl b/test/jump.jl index 8b730a1..1207683 100644 --- a/test/jump.jl +++ b/test/jump.jl @@ -1,4 +1,35 @@ import JuMP +import QUBODrivers +import QUBOTools + +struct FakePEPSTopology <: TenSolver.AbstractStructuredTopology + variables::Int +end + +function TenSolver._solve_ising(backend::TenSolver.PEPSBackend{FakePEPSTopology}, model::TenSolver.IsingModel; kwargs...) + @test backend.topology.variables == length(model.h) + @test backend.beta == 1.75 + @test backend.bond_dim == 3 + @test backend.max_states == 4 + @test backend.cutoff_prob == 0.0 + @test !backend.onGPU + @test backend.contraction == :svd + @test backend.num_sweeps == 1 + @test backend.transformations == :identity + @test backend.local_dimension == 2 + + states = [[1, 1], [1, 0]] + energies = [-3.0, -1.0] + probabilities = [0.8, 0.2] + metadata = Dict{String, Any}( + "backend" => "FakePEPS", + "topology" => "fake", + "selected_transformation" => "identity", + "largest_discarded_probability" => 0.0, + ) + + return first(energies), TenSolver.PEPSSolution{Float64}(states, energies, probabilities, metadata, nothing) +end @testset "JuMP interface" begin dim = 5 @@ -22,3 +53,219 @@ import JuMP # Same solution @test JuMP.value.(x) == x0 end + +@testset "JuMP backend attributes" begin + @testset "Default backend" begin + model = JuMP.Model(TenSolver.Optimizer) + JuMP.set_silent(model) + @JuMP.variable(model, x, Bin) + @JuMP.objective(model, Min, -x) + + JuMP.optimize!(model) + + @test JuMP.objective_value(model) ≈ -1.0 + @test JuMP.value(x) ≈ 1.0 + end + + @testset "Explicit DMRG backend" begin + model = JuMP.Model(TenSolver.Optimizer) + JuMP.set_silent(model) + JuMP.set_attribute(model, "backend", :dmrg) + JuMP.set_attribute(model, "num_reads", 4) + @JuMP.variable(model, x, Bin) + @JuMP.objective(model, Min, -x) + + JuMP.optimize!(model) + + @test JuMP.objective_value(model) ≈ -1.0 + @test JuMP.value(x) ≈ 1.0 + end + + @testset "String backend normalization" begin + model = JuMP.Model(TenSolver.Optimizer) + JuMP.set_silent(model) + JuMP.set_attribute(model, "backend", "dmrg") + @JuMP.variable(model, x, Bin) + @JuMP.objective(model, Min, -x) + + JuMP.optimize!(model) + + @test JuMP.objective_value(model) ≈ -1.0 + end + + @testset "PEPS requires topology metadata" begin + model = JuMP.Model(TenSolver.Optimizer) + JuMP.set_silent(model) + JuMP.set_attribute(model, "backend", :peps) + @JuMP.variable(model, x, Bin) + @JuMP.objective(model, Min, -x) + + err = try + JuMP.optimize!(model) + catch err + err + end + + @test err isa ArgumentError + @test occursin("peps_topology", sprint(showerror, err)) + @test occursin("(m, n)", sprint(showerror, err)) + end + + @testset "Unavailable PEPS extension errors clearly" begin + has_spinglasspeps_components = all(pkg -> !isnothing(Base.find_package(pkg)), ( + "SpinGlassNetworks", + "SpinGlassEngine", + "SpinGlassTensors", + )) + + if has_spinglasspeps_components + @test_skip "SpinGlassPEPS component packages are available; unavailable-extension error path does not apply." + else + model = JuMP.Model(TenSolver.Optimizer) + JuMP.set_silent(model) + JuMP.set_attribute(model, "backend", :peps) + JuMP.set_attribute(model, "peps_layout", :square) + JuMP.set_attribute(model, "peps_topology", (1, 1)) + JuMP.set_attribute(model, "peps_beta", 1.5) + JuMP.set_attribute(model, "peps_bond_dim", 4) + JuMP.set_attribute(model, "peps_max_states", 2) + JuMP.set_attribute(model, "peps_strategy", :svd) + @JuMP.variable(model, x, Bin) + @JuMP.objective(model, Min, -x) + + err = try + JuMP.optimize!(model) + catch err + err + end + + @test err isa ArgumentError + @test occursin("PEPSBackend is not available", sprint(showerror, err)) + @test occursin("SpinGlassNetworks", sprint(showerror, err)) + end + end + + @testset "Optional PEPS optimizer solve" begin + has_spinglasspeps_components = all(pkg -> !isnothing(Base.find_package(pkg)), ( + "SpinGlassNetworks", + "SpinGlassEngine", + "SpinGlassTensors", + )) + + if !has_spinglasspeps_components + @test_skip "SpinGlassPEPS component packages are not available in this environment." + else + import SpinGlassNetworks + import SpinGlassEngine + import SpinGlassTensors + + Q = [ + -1.0 0.5 0.0 0.0 + 0.0 -0.5 0.0 0.0 + 0.0 0.0 -0.25 0.25 + 0.0 0.0 0.0 -0.75 + ] + l = [0.0, 0.25, -0.25, 0.0] + c = 0.125 + objective(x) = dot(x, Q, x) + dot(l, x) + c + exact_energy, _ = brute_force(objective, Float64, 4) + + model = JuMP.Model(TenSolver.Optimizer) + JuMP.set_silent(model) + JuMP.set_attribute(model, "backend", :peps) + JuMP.set_attribute(model, "peps_layout", :square) + JuMP.set_attribute(model, "peps_topology", (2, 2)) + JuMP.set_attribute(model, "peps_beta", 2.0) + JuMP.set_attribute(model, "peps_bond_dim", 4) + JuMP.set_attribute(model, "peps_max_states", 4) + JuMP.set_attribute(model, "peps_cutoff_prob", 0.0) + JuMP.set_attribute(model, "peps_strategy", :svd) + JuMP.set_attribute(model, "peps_transformations", :identity) + @JuMP.variable(model, x[1:4], Bin) + @JuMP.objective( + model, + Min, + sum(Q[i, j] * x[i] * x[j] for i in 1:4, j in 1:4) + + sum(l[i] * x[i] for i in 1:4) + c + ) + + JuMP.optimize!(model) + + state = round.(Int, JuMP.value.(x)) + @test JuMP.objective_value(model) ≈ exact_energy atol = 1e-6 + @test objective(state) ≈ JuMP.objective_value(model) atol = 1e-6 + end + end + + @testset "Fake PEPS optimizer solve" begin + model = JuMP.Model(TenSolver.Optimizer) + JuMP.set_silent(model) + JuMP.set_attribute(model, "backend", :peps) + JuMP.set_attribute(model, "peps_topology", FakePEPSTopology(2)) + JuMP.set_attribute(model, "peps_beta", 1.75) + JuMP.set_attribute(model, "peps_bond_dim", 3) + JuMP.set_attribute(model, "peps_max_states", 4) + JuMP.set_attribute(model, "peps_cutoff_prob", 0.0) + JuMP.set_attribute(model, "peps_strategy", :svd) + JuMP.set_attribute(model, "peps_transformations", :identity) + JuMP.set_attribute(model, "peps_local_dimension", 2) + JuMP.set_attribute(model, QUBODrivers.FinalNumberOfReads(), 5) + @JuMP.variable(model, x[1:2], Bin) + @JuMP.objective(model, Min, -x[1] - 2x[2]) + + JuMP.optimize!(model) + + solution = QUBOTools.solution(JuMP.unsafe_backend(model)) + metadata = QUBOTools.metadata(solution) + + @test JuMP.objective_value(model) ≈ -3.0 + @test round.(Int, JuMP.value.(x)) == [1, 1] + @test QUBOTools.reads(solution) == 5 + @test QUBOTools.state(solution, 1) == [1, 1] + @test QUBOTools.reads(solution, 1) == 4 + @test QUBOTools.state(solution, 2) == [1, 0] + @test QUBOTools.reads(solution, 2) == 1 + @test isempty(QUBODrivers.validate_metadata(solution)) + @test metadata["algorithm"]["name"] == "FakePEPS" + @test metadata["backend"]["name"] == "TenSolver" + @test metadata["backend"]["version"] == TenSolver.__VERSION__ + @test metadata["optimizer"]["evaluations"] == 2 + @test metadata["reads"]["final_number_of_reads"] == 5 + @test metadata["tensolver"]["peps"]["topology"] == "fake" + @test metadata["tensolver"]["peps"]["candidate_states"] == 2 + end + + @testset "PEPS SampleSet adaptation" begin + peps = TenSolver.PEPSSolution{Float64}( + [[1, 0], [0, 1]], + [-2.0, -1.0], + [0.75, 0.25], + Dict{String, Any}( + "backend" => "SpinGlassPEPS", + "topology" => "square", + "selected_transformation" => "identity", + "largest_discarded_probability" => 0.01, + ), + nothing, + ) + Q = [0.0 -1.0; 0.0 0.0] + l = [0.0, -0.5] + samples = TenSolver._qubo_samples(Float64, peps, l, Q, 1.0, 0.0, 3) + + @test getfield.(samples, :state) == [[1, 0], [0, 1]] + @test getfield.(samples, :value) == [0.0, -0.5] + @test getfield.(samples, :reads) == [2, 1] + + metadata = Dict{String, Any}( + "origin" => "TenSolver", + "time" => Dict{String, Any}("effective" => 1.25), + ) + TenSolver._add_backend_metadata!(metadata, peps) + + @test metadata["tensolver"]["peps"]["backend"] == "SpinGlassPEPS" + @test metadata["tensolver"]["peps"]["topology"] == "square" + @test metadata["tensolver"]["peps"]["candidate_states"] == 2 + @test metadata["tensolver"]["peps"]["effective_time"] == 1.25 + @test metadata["tensolver"]["peps"]["largest_discarded_probability"] == 0.01 + end +end