diff --git a/Project.toml b/Project.toml index 58be85a..d9cbf8d 100644 --- a/Project.toml +++ b/Project.toml @@ -2,23 +2,25 @@ name = "QuantumStateTransfer" uuid = "5afda8b7-781f-4fb2-840f-bdb502253c46" keywords = ["quantum computing", "quantum information", "linear algebra", "graph theory", "optimization"] license = "MIT" +version = "0.1.0-dev" authors = ["Luis M. B. Varona ", "Nathaniel Johnston "] description = "A Julia toolbox for modeling state transfer on quantum networks." homepage = "https://graphquantum.github.io/QuantumStateTransfer.jl/" maintainers = ["Luis M. B. Varona "] readme = "README.md" repository = "https://github.com/GraphQuantum/QuantumStateTransfer.jl" -version = "0.1.0-dev" [deps] DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" [compat] DataStructures = "0.18.15 - 0.19" Graphs = "1.10" LinearAlgebra = "1.10" +Optim = "1.13.2" PrecompileTools = "1.2" julia = "1.10" diff --git a/src/EpsilonOptimization/EpsilonOptimization.jl b/src/EpsilonOptimization/EpsilonOptimization.jl new file mode 100644 index 0000000..a880a00 --- /dev/null +++ b/src/EpsilonOptimization/EpsilonOptimization.jl @@ -0,0 +1,39 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +""" + EpsilonOptimization + +Optimization algorithms with guaranteed finite ``ϵ``-convergence to the true global minimum +of ``ℝⁿ → ℝ`` functions over hyperrectangular domains. +""" +module EpsilonOptimization + +using QuantumStateTransfer: NotImplementedError + +using DataStructures: BinaryMinHeap +using LinearAlgebra: norm +using Optim + +export + # Types + AbstractEpsilonSolver, + EpsilonMinimizationResult, + + # Core functions + epsilon_minimize, + + # Solvers + LipschitzBranchAndBound, + AlphaBranchAndBound + +include("types.jl") +include("core.jl") + +include("solvers/lipschitz_branch_and_bound.jl") +include("solvers/alpha_branch_and_bound.jl") + +end diff --git a/src/EpsilonOptimization/core.jl b/src/EpsilonOptimization/core.jl new file mode 100644 index 0000000..e2285e6 --- /dev/null +++ b/src/EpsilonOptimization/core.jl @@ -0,0 +1,46 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +""" + epsilon_minimize(f, lower, upper, solver) + +Converge to within an arbitrary ``ϵ`` of the true global minimum of an ``ℝⁿ → ℝ`` function +`f` over the hyperrectangular domain defined by the bounds `lower` and `upper`, using the +specified ``ϵ``-optimization `solver`. + +# Arguments +- `f::Function`: The objective function to be minimized. Must map from `ℝⁿ → ℝ`. +- `lower::Tx<:AbstractVector{<:Real}`: The lower bounds of the hyperrectangular domain. +- `upper::Tx<:AbstractVector{<:Real}`: The upper bounds of the hyperrectangular domain. +- `solver::AbstractEpsilonSolver`: The ``ϵ``-optimization solver to use. + +# Returns +- `::EpsilonMinimizationResult{Tx,<:Real}}`: The result of the ``ϵ``-optimization. + +# Examples +[TODO: Write here] +""" +function epsilon_minimize( + f::Function, lower::Tx, upper::Tx, solver::S +) where {Tx<:AbstractVector{<:Real},S<:AbstractEpsilonSolver} + if length(lower) != length(upper) + throw(ArgumentError("Lower and upper bounds must have the same dimension")) + end + + if any(lower .> upper) + throw( + ArgumentError("Lower bound must be entrywise less than or equal to upper bound") + ) + end + + return _epsilon_minimize_impl(f, lower, upper, solver) +end + +function _epsilon_minimize_impl( + ::Function, ::Tx, ::Tx, ::S +) where {Tx<:AbstractVector{<:Real},S<:AbstractEpsilonSolver} + throw(NotImplementedError(_epsilon_minimize_impl, :solver, S, AbstractEpsilonSolver)) +end diff --git a/src/EpsilonOptimization/solvers/alpha_branch_and_bound.jl b/src/EpsilonOptimization/solvers/alpha_branch_and_bound.jl new file mode 100644 index 0000000..c003c32 --- /dev/null +++ b/src/EpsilonOptimization/solvers/alpha_branch_and_bound.jl @@ -0,0 +1,147 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +""" + AlphaBranchAndBound <: AbstractEpsilonSolver + +[TODO: Write here] +""" +struct AlphaBranchAndBound <: AbstractEpsilonSolver + epsilon::Real + threshold::Union{<:Real,Nothing} + max_iterations::Union{<:Integer,Nothing} + alpha::Real +end + +""" + ABBHyperrectangle{Tx,Tf} + +[TODO: Write here] +""" +struct ABBHyperrectangle{Tx<:AbstractVector{<:Real},Tf<:Real} + lower::Tx + upper::Tx + x_min::Tx + lower_bound::Tf + + function ABBHyperrectangle( + lower::Tx, upper::Tx, f::Function, alpha::Real + ) where {Tx<:AbstractVector{<:Real}} + f_convex(x::Tx) = f(x) + alpha * sum((x .- lower) .* (upper .- x)) + x0 = lower .+ (upper .- lower) ./ 2 + + res = optimize(f_convex, lower, upper, x0, Fminbox(LBFGS())) + x_min = res.minimizer + lower_bound = res.minimum + + return new{Tx,typeof(lower_bound)}(lower, upper, x_min, lower_bound) + end +end + +function Base.isless(rect1::ABBHyperrectangle, rect2::ABBHyperrectangle) + return rect1.lower_bound < rect2.lower_bound +end + +function _epsilon_minimize_impl( + f::Function, lower::Tx, upper::Tx, solver::AlphaBranchAndBound +) where {Tx<:AbstractVector{<:Real}} + epsilon = solver.epsilon + alpha = solver.alpha + + if isnothing(solver.threshold) + threshold = -Inf + else + threshold = solver.threshold + end + + if isnothing(solver.max_iterations) + max_iterations = Inf + else + max_iterations = solver.max_iterations + end + + rect_init = ABBHyperrectangle(lower, upper, f, alpha) + rects_cand = BinaryMinHeap{ABBHyperrectangle{Tx,typeof(rect_init.lower_bound)}}() + push!(rects_cand, rect_init) + + minimizer = rect_init.x_min + minimum = f(rect_init.x_min) + lower_bound = Inf + iterations = 0 + + while ( + iterations < max_iterations && + !isempty(rects_cand) && + min(lower_bound, first(rects_cand).lower_bound) <= threshold && + minimum - threshold >= epsilon + ) + rect = pop!(rects_cand) + f_val = f(rect.x_min) + + if f_val < minimum + minimizer = rect.x_min + minimum = f_val + end + + children = _abb_split_hyperrectangle(rect, f, alpha) + + for child in children + f_val_child = f(child.x_min) + + if f_val_child < minimum + minimizer = child.x_min + minimum = f_val_child + end + + if minimum - child.lower_bound >= epsilon + push!(rects_cand, child) + else + lower_bound = min(lower_bound, child.lower_bound) + end + end + + iterations += 1 + end + + if !isempty(rects_cand) + lower_bound = min(lower_bound, first(rects_cand).lower_bound) + end + + return EpsilonMinimizationResult( + solver, + epsilon, + lower, + upper, + minimizer, + minimum, + ( + epsilon_optimal=minimum - lower_bound < epsilon, + threshold_reached=minimum - threshold < epsilon, + threshold_unreachable=lower_bound > threshold, + iterations=iterations >= max_iterations, + ), + ) +end + +function _abb_split_hyperrectangle(rect::ABBHyperrectangle, f::Function, alpha::Real) + lower = rect.lower + upper = rect.upper + + dim_split = argmax(upper .- lower) + mid = lower[dim_split] + (upper[dim_split] - lower[dim_split]) / 2 + + lower1 = copy(lower) + upper1 = copy(upper) + upper1[dim_split] = mid + child1 = ABBHyperrectangle(lower1, upper1, f, alpha) + + lower2 = copy(lower) + upper2 = copy(upper) + lower2[dim_split] = mid + child2 = ABBHyperrectangle(lower2, upper2, f, alpha) + + return child1, child2 +end diff --git a/src/EpsilonOptimization/solvers/lipschitz_branch_and_bound.jl b/src/EpsilonOptimization/solvers/lipschitz_branch_and_bound.jl new file mode 100644 index 0000000..df52346 --- /dev/null +++ b/src/EpsilonOptimization/solvers/lipschitz_branch_and_bound.jl @@ -0,0 +1,144 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +""" + LipschitzBranchAndBound <: AbstractEpsilonSolver + +[TODO: Write here] +""" +struct LipschitzBranchAndBound <: AbstractEpsilonSolver + epsilon::Real + threshold::Union{<:Real,Nothing} + max_iterations::Union{<:Integer,Nothing} + lipschitz_constant::Real +end + +""" + LBBHyperrectangle{Tx,Tf} + +[TODO: Write here] +""" +struct LBBHyperrectangle{Tx<:AbstractVector{<:Real},Tf<:Real} + lower::Tx + upper::Tx + center::Tx + f_center::Tf + lower_bound::Tf + + function LBBHyperrectangle( + lower::Tx, upper::Tx, f::Function, lipschitz_constant::Real + ) where {Tx<:AbstractVector{<:Real}} + center = lower .+ (upper .- lower) ./ 2 + f_center = f(center) + diameter = norm(upper .- lower) + lower_bound = f_center - (lipschitz_constant / 2) * diameter + return new{Tx,typeof(f_center)}(lower, upper, center, f_center, lower_bound) + end +end + +function Base.isless(rect1::LBBHyperrectangle, rect2::LBBHyperrectangle) + return rect1.lower_bound < rect2.lower_bound +end + +function _epsilon_minimize_impl( + f::Function, lower::Tx, upper::Tx, solver::LipschitzBranchAndBound +) where {Tx<:AbstractVector{<:Real}} + epsilon = solver.epsilon + lipschitz_constant = solver.lipschitz_constant + + if isnothing(solver.threshold) + threshold = -Inf + else + threshold = solver.threshold + end + + if isnothing(solver.max_iterations) + max_iterations = Inf + else + max_iterations = solver.max_iterations + end + + rect_init = LBBHyperrectangle(lower, upper, f, solver.lipschitz_constant) + rects_cand = BinaryMinHeap{LBBHyperrectangle{Tx,typeof(rect_init.f_center)}}() + push!(rects_cand, rect_init) + + minimizer = rect_init.center + minimum = rect_init.f_center + lower_bound = Inf + iterations = 0 + + while ( + iterations < max_iterations && + !isempty(rects_cand) && + min(lower_bound, first(rects_cand).lower_bound) <= threshold && + minimum - threshold >= epsilon + ) + rect = pop!(rects_cand) + + if rect.f_center < minimum + minimizer = rect.center + minimum = rect.f_center + end + + children = _lbb_split_hyperrectangle(rect, f, lipschitz_constant) + + for child in children + if child.f_center < minimum + minimizer = child.center + minimum = child.f_center + end + + if minimum - child.lower_bound >= epsilon + push!(rects_cand, child) + else + lower_bound = min(lower_bound, child.lower_bound) + end + end + + iterations += 1 + end + + if !isempty(rects_cand) + lower_bound = min(lower_bound, first(rects_cand).lower_bound) + end + + return EpsilonMinimizationResult( + solver, + epsilon, + lower, + upper, + minimizer, + minimum, + ( + epsilon_optimal=minimum - lower_bound < epsilon, + threshold_reached=minimum - threshold < epsilon, + threshold_unreachable=lower_bound > threshold, + iterations=iterations >= max_iterations, + ), + ) +end + +function _lbb_split_hyperrectangle( + rect::LBBHyperrectangle, f::Function, lipschitz_constant::Real +) + lower = rect.lower + upper = rect.upper + + dim_split = argmax(upper .- lower) + mid = lower[dim_split] + (upper[dim_split] - lower[dim_split]) / 2 + + lower1 = copy(lower) + upper1 = copy(upper) + upper1[dim_split] = mid + child1 = LBBHyperrectangle(lower1, upper1, f, lipschitz_constant) + + lower2 = copy(lower) + upper2 = copy(upper) + lower2[dim_split] = mid + child2 = LBBHyperrectangle(lower2, upper2, f, lipschitz_constant) + + return child1, child2 +end diff --git a/src/EpsilonOptimization/types.jl b/src/EpsilonOptimization/types.jl new file mode 100644 index 0000000..5de5ad7 --- /dev/null +++ b/src/EpsilonOptimization/types.jl @@ -0,0 +1,80 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +""" + AbstractEpsilonSolver + +Abstract base type for all ``ϵ``-optimization solvers. + +# Interface +Concrete subtypes of `AbstractEpsilonSolver` must implement the following methods: +- `Base.summary(::T) where {T<:AbstractAlgorithm}`: returns a `String` indicating the name + of the algorithm (e.g., `"Lipschitz branch-and-bound"`). + +Concrete subtypes of `AbstractEpsilonSolver` must also have the following fields: +- `epsilon::Real`: The tolerance for ``ϵ``-convergence. +- `threshold::Union{<:Real,Nothing}`: An optional threshold value resulting in termination + once either (1) a function value less than `threshold + epsilon` is found or (2) it is + determined that no function value less than or equal to `threshold` exists. +- `max_iterations::Union{<:Integer,Nothing}`: An optional maximum number of iterations after + which the algorithm will terminate. +""" +abstract type AbstractEpsilonSolver end + +""" + EpsilonMinimizationResult{Tx,Tf} + +Output struct for ``ϵ``-optimization results. + +# Fields +- `algorithm::AbstractEpsilonSolver`: The algorithm used for the optimization. +- `lower::Tx<:AbstractVector{<:Real}}`: The lower bounds of the hyperrectangular domain. +- `upper::Tx<:AbstractVector{<:Real}}`: The upper bounds of the hyperrectangular domain. +- `epsilon::Real`: The tolerance provided for ``ϵ``-convergence. +- `minimizer::Tx<:AbstractVector{<:Real}}`: The estimated location of the global minimizer. +- `minimum::Tf<:Real`: The estimated global minimum function value. +- `stopped_by::NamedTuple`: The reason(s) for termination, namely: + - `:epsilon_optimal::Bool`: Whether ``ϵ``-optimality was achieved. + - `:threshold_reached::Bool`: Whether ta function value less than `threshold + epsilon` + was found for some optionally provided `threshold`. + - `:threshold_unreachable::Bool`: Whether it was determined that no function value less + than or equal to `threshold` existed for some optionally provided `threshold`. + - `:iterations::Bool`: Whether termination occurred due to reaching the maximum number + of iterations for some optionally provided `max_iterations`. +""" +struct EpsilonMinimizationResult{Tx<:AbstractVector{<:Real},Tf<:Real} + algorithm::AbstractEpsilonSolver + lower::Tx + upper::Tx + epsilon::Real + minimizer::Tx + minimum::Tf + stopped_by::NamedTuple{ + (:epsilon_optimal, :threshold_reached, :threshold_unreachable, :iterations), + Tuple{Bool,Bool,Bool,Bool}, + } +end + +# The `Base.show` override here takes heavy inspiration from the Optim.jl package +function Base.show(io::IO, res::EpsilonMinimizationResult) + println(io, "Results of Epsilon Minimization Algorithm") + print(io, " * Status: ") + + if res.stopped_by.epsilon_optimal || res.stopped_by.threshold_reached + println(io, "success") + else + println(io, "failure") + end + + println(io, " * Algorithm: $(summary(res.algorithm))") + println(io, " * Lower bounds: $(res.lower)") + println(io, " * Upper bounds: $(res.upper)") + println(io, " * Epsilon tolerance: $(res.epsilon)") + println(io, " * Estimated minimizer: $(res.minimizer)") + println(io, " * Estimated minimum: $(res.minimum)") + + return nothing +end diff --git a/src/QuantumStateTransfer.jl b/src/QuantumStateTransfer.jl index fb37123..8b36be3 100644 --- a/src/QuantumStateTransfer.jl +++ b/src/QuantumStateTransfer.jl @@ -4,8 +4,28 @@ # http://opensource.org/licenses/MIT>. This file may not be copied, modified, or # distributed except according to those terms. +""" + QuantumStateTransfer + +A Julia toolbox for modeling state transfer on quantum networks. +""" module QuantumStateTransfer -# Write your package code here. +using DataStructures +using Graphs +using LinearAlgebra +using PrecompileTools: @setup_workload, @compile_workload + +include("utils.jl") +include("types.jl") + +include("EpsilonOptimization/EpsilonOptimization.jl") +using .EpsilonOptimization + +include("core.jl") + +# TODO: Exports + +include("startup.jl") end diff --git a/src/core.jl b/src/core.jl new file mode 100644 index 0000000..fa351df --- /dev/null +++ b/src/core.jl @@ -0,0 +1,7 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +# TODO: Write here diff --git a/src/startup.jl b/src/startup.jl new file mode 100644 index 0000000..fa351df --- /dev/null +++ b/src/startup.jl @@ -0,0 +1,7 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +# TODO: Write here diff --git a/src/types.jl b/src/types.jl new file mode 100644 index 0000000..018e7fc --- /dev/null +++ b/src/types.jl @@ -0,0 +1,86 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +""" + NotImplementedError{Nothing}(f, subtype, abstracttype) + NotImplementedError{Symbol}(f, arg, subtype, abstracttype) + +An exception indicating that a function lacks dispatch to handle a specific argument type. + +Semantically, this differs from `MethodError` in that it connotes a developer-side failure +to implement a method rather than erroneous user input. Throughout this package, it is often +used to warn when an existing function with multiple dispatch on some abstract type is +called on a newly created subtype for which no method has been defined. + +# Fields +- `f::Function`: the function called. +- `arg::Symbol`: the name of the argument with the unsupported type, if the function has + multiple arguments. If the function has only one argument, this field should be set to + `nothing`. +- `subtype::Type`: the type of the argument. May be the actual concrete type or some + intermediate supertype. (For instance, if the relevant input has concrete type `A` with + hierarchy `A <: B <: C` and the `abstracttype` field is `C`, then both `A` and `B` are + perfectly valid choices for `subtype`.) +- `abstracttype::Type`: the abstract type under which the argument is meant to fall. + +# Constructors +- `NotImplementedError(::Function, ::Type, ::Type)`: constructs a new `NotImplementedError` + instance for a single-argument function. Throws an error if the second type is not + abstract or the first type is not a subtype of the second. +- `NotImplementedError(::Function, ::Symbol, ::Type, ::Type)`: constructs a new + `NotImplementedError` instance for a multi-argument function. Throws an error if the + second type is not abstract or the first type is not a subtype of the second. + +# Supertype Hierarchy +`NotImplementedError` <: `Exception` + +# Notes +This struct was taken from one of the authors' other packages, MatrixBandwidth.jl. +""" +struct NotImplementedError{T<:Union{Nothing,Symbol}} <: Exception + f::Function + arg::T + subtype::Type + abstracttype::Type + + function NotImplementedError(f::Function, subtype::Type, abstracttype::Type) + return NotImplementedError(f, nothing, subtype, abstracttype) + end + + function NotImplementedError( + f::Function, arg::T, subtype::Type, abstracttype::Type + ) where {T<:Union{Nothing,Symbol}} + if !isabstracttype(abstracttype) + throw(ArgumentError("Expected an abstract type, got $abstracttype")) + end + + if !(subtype <: abstracttype) + throw(ArgumentError("Expected a subtype of $abstracttype, got $subtype")) + end + + return new{T}(f, arg, subtype, abstracttype) + end +end + +function Base.showerror(io::IO, e::NotImplementedError{Nothing}) + print( + io, + """NotImplementedError with $(e.subtype): + $(e.f) is not yet implemented for this subtype of $(e.abstracttype). + Try defining method dispatch manually if this is a newly created subtype.""", + ) + return nothing +end + +function Base.showerror(io::IO, e::NotImplementedError{Symbol}) + print( + io, + """NotImplementedError with argument $(e.arg)::$(e.subtype): + $(e.f) is not yet implemented for this subtype of $(e.abstracttype). + Try defining method dispatch manually if this is a newly created subtype.""", + ) + return nothing +end diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..fa351df --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,7 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +# TODO: Write here