From 72751b988f337cfb43832942b79db66700d60a4a Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Thu, 21 Dec 2023 11:07:35 +0100 Subject: [PATCH 1/2] make the bounds type open to extension --- Project.toml | 2 +- docs/src/metabolic-modeling.jl | 14 ++++---- docs/src/quadratic-optimization.jl | 20 +++++------ src/bound.jl | 55 ++++++++++++++++++++++++++---- src/constraint.jl | 42 +++++++++++------------ src/constraint_tree.jl | 8 ++--- test/misc.jl | 8 +++-- 7 files changed, 94 insertions(+), 55 deletions(-) diff --git a/Project.toml b/Project.toml index 23a3312..4bae20b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ConstraintTrees" uuid = "5515826b-29c3-47a5-8849-8513ac836620" authors = ["The developers of ConstraintTrees.jl"] -version = "0.6.1" +version = "0.7.0" [deps] DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" diff --git a/docs/src/metabolic-modeling.jl b/docs/src/metabolic-modeling.jl index 0feaaa9..e1f9aa4 100644 --- a/docs/src/metabolic-modeling.jl +++ b/docs/src/metabolic-modeling.jl @@ -226,12 +226,12 @@ function optimized_vars(cs::C.ConstraintTree, objective::C.LinearValue, optimize JuMP.@objective(model, JuMP.MAX_SENSE, C.substitute(objective, x)) function add_constraint(c::C.Constraint) b = c.bound - if b isa Float64 - JuMP.@constraint(model, C.substitute(c.value, x) == b) - elseif b isa Tuple{Float64,Float64} + if b isa C.EqualTo + JuMP.@constraint(model, C.substitute(c.value, x) == b.equal_to) + elseif b isa C.Between val = C.substitute(c.value, x) - isinf(b[1]) || JuMP.@constraint(model, val >= b[1]) - isinf(b[2]) || JuMP.@constraint(model, val <= b[2]) + isinf(b.lower) || JuMP.@constraint(model, val >= b.lower) + isinf(b.upper) || JuMP.@constraint(model, val <= b.upper) end end function add_constraint(c::C.ConstraintTree) @@ -346,9 +346,9 @@ Dict(k => v.fluxes.R_BIOMASS_Ecoli_core_w_GAM for (k, v) in result.community) # both dot-access and array-index syntax. # You can thus, e.g., set a single bound: -c.exchanges.oxygen.bound = (-20.0, 20.0) +c.exchanges.oxygen.bound = C.Between(-20.0, 20.0) -# ...or rebuild a whole constraint: +# ...or rebuild a whole constraint (using a tuple shortcut for [`Between`](@ref)): c.exchanges.biomass = C.Constraint(c.exchanges.biomass.value, (-20, 20)) # ...or even add new constraints, here using the index syntax for demonstration: diff --git a/docs/src/quadratic-optimization.jl b/docs/src/quadratic-optimization.jl index 9728228..b976fca 100644 --- a/docs/src/quadratic-optimization.jl +++ b/docs/src/quadratic-optimization.jl @@ -32,7 +32,7 @@ error_val = # This allows us to naturally express quadratic constraint (e.g., that an error # must not be too big); and directly observe the error values in the system. -system = :vars^system * :error^C.Constraint(error_val, (0, 100)) +system = :vars^system * :error^C.Constraint(error_val, C.Between(0, 100)) # (For simplicity, you can also use the `Constraint` constructor to make # quadratic constraints out of `QuadraticValue`s -- it will overload properly.) @@ -58,7 +58,7 @@ ellipse_system = C.ConstraintTree( :point => point, :in_area => C.Constraint( C.squared(point.x.value) / 4 + C.squared(10.0 - point.y.value), - (-Inf, 1.0), + C.Between(-Inf, 1.0), ), ) @@ -95,22 +95,18 @@ s *= # translates the constraints into JuMP `Model`s to support the quadratic # constraints. import JuMP -function optimized_vars( - cs::C.ConstraintTree, - objective::Union{C.LinearValue,C.QuadraticValue}, - optimizer, -) +function optimized_vars(cs::C.ConstraintTree, objective::C.Value, optimizer) model = JuMP.Model(optimizer) JuMP.@variable(model, x[1:C.var_count(cs)]) JuMP.@objective(model, JuMP.MAX_SENSE, C.substitute(objective, x)) function add_constraint(c::C.Constraint) b = c.bound - if b isa Float64 - JuMP.@constraint(model, C.substitute(c.value, x) == b) - elseif b isa Tuple{Float64,Float64} + if b isa C.EqualTo + JuMP.@constraint(model, C.substitute(c.value, x) == b.equal_to) + elseif b isa C.Between val = C.substitute(c.value, x) - isinf(b[1]) || JuMP.@constraint(model, val >= b[1]) - isinf(b[2]) || JuMP.@constraint(model, val <= b[2]) + isinf(b.lower) || JuMP.@constraint(model, val >= b.lower) + isinf(b.upper) || JuMP.@constraint(model, val <= b.upper) end end function add_constraint(c::C.ConstraintTree) diff --git a/src/bound.jl b/src/bound.jl index 3bce0dc..10a4c92 100644 --- a/src/bound.jl +++ b/src/bound.jl @@ -2,16 +2,59 @@ """ $(TYPEDEF) -Convenience shortcut for "interval" bound; consisting of lower and upper bound +Abstract type of all bounds usable in constraints, including [`Between`](@ref) +and [`EqualTo`](@ref). +""" +abstract type Bound end + +""" +$(TYPEDEF) + +Representation of an "equality" bound; contains the single "equal to this" value. + +# Fields +$(TYPEDFIELDS) """ -const IntervalBound = Tuple{Float64,Float64} +Base.@kwdef mutable struct EqualTo <: Bound + "Equality bound value" + equal_to::Float64 + + EqualTo(x::Real) = new(Float64(x)) +end + +Base.:-(x::EqualTo) = -1 * x +Base.:*(a::EqualTo, b::Real) = b * a +Base.:/(a::EqualTo, b::Real) = EqualTo(a.equal_to / b) +Base.:*(a::Real, b::EqualTo) = EqualTo(a * b.equal_to) + +""" +$(TYPEDEF) + +Representation of an "interval" bound; consisting of lower and upper bound +value. + +# Fields +$(TYPEDFIELDS) +""" +Base.@kwdef mutable struct Between <: Bound + "Lower bound" + lower::Float64 = -Inf + "Upper bound" + upper::Float64 = Inf + + Between(x::Real, y::Real) = new(Float64(x), Float64(y)) +end + +Base.:-(x::Between) = -1 * x +Base.:*(a::Between, b::Real) = b * a +Base.:/(a::Between, b::Real) = Between(a.lower / b, a.upper / b) +Base.:*(a::Real, b::Between) = Between(a * b.lower, a * b.upper) """ $(TYPEDEF) -Shortcut for possible bounds: either no bound is present (`nothing`), or a -single number is interpreted as an exact equality bound, or a tuple of 2 -numbers is interpreted as an interval bound. +Shortcut for all possible [`Bound`](@ref)s including the "empty" bound that +does not constraint anything (represented by `nothing`). """ -const Bound = Union{Nothing,Float64,IntervalBound} +const MaybeBound = Union{Nothing,Bound} diff --git a/src/constraint.jl b/src/constraint.jl index 2459fab..7c8d941 100644 --- a/src/constraint.jl +++ b/src/constraint.jl @@ -2,10 +2,8 @@ """ $(TYPEDEF) -A representation of a single constraint that limits a given value by a specific -[`Bound`](@ref). - -Constraints may be scaled linearly, i.e., multiplied by real-number constants. +A representation of a single constraint that may limit the given value by a +specific [`Bound`](@ref). Constraints without a bound (`nothing` in the `bound` field) are possible; these have no impact on the optimization problem but the associated `value` @@ -14,34 +12,32 @@ becomes easily accessible for inspection and building other constraints. # Fields $(TYPEDFIELDS) """ -Base.@kwdef mutable struct Constraint{V} +Base.@kwdef mutable struct Constraint{V,B} "A value (typically a [`LinearValue`](@ref) or a [`QuadraticValue`](@ref)) that describes what the constraint constraints." value::V - "A bound that the `value` must satisfy." - bound::Bound = nothing + "A bound that the `value` must satisfy. Should be a subtype of + [`MaybeBound`](@ref): Either `nothing` if there's no bound, or e.g. + [`EqualTo`](@ref), [`Between`](@ref) or similar structs." + bound::B = nothing - function Constraint(v::T, b::Bound = nothing) where {T<:Value} - new{T}(v, b) + function Constraint(v::T, b::U = nothing) where {T<:Value,U<:MaybeBound} + new{T,U}(v, b) end end -Constraint(v::T, b::Int) where {T<:Value} = Constraint(v, Float64(b)) +Constraint(v::T, b::Real) where {T<:Value} = Constraint(v, EqualTo(b)) Constraint(v::T, b::Tuple{X,Y}) where {T<:Value,X<:Real,Y<:Real} = - Constraint(v, Float64.(b)) + Constraint(v, Between(Float64.(b)...)) -Base.:-(a::Constraint) = -1 * a -Base.:*(a::Real, b::Constraint) = b * a -Base.:*(a::Constraint, b::Real) = Constraint( - value = a.value * b, - bound = a.bound isa Float64 ? a.bound * b : - a.bound isa Tuple{Float64,Float64} ? a.bound .* b : nothing, -) -Base.:/(a::Constraint, b::Real) = Constraint( - value = a.value / b, - bound = a.bound isa Float64 ? a.bound / b : - a.bound isa Tuple{Float64,Float64} ? a.bound ./ b : nothing, -) +Base.:-(a::Constraint) = + Constraint(value = -a.value, bound = isnothing(a.bound) ? nothing : -a.bound) +Base.:*(a::Real, b::Constraint) = + Constraint(value = a * b.value, bound = isnothing(b.bound) ? nothing : a * b.bound) +Base.:*(a::Constraint, b::Real) = + Constraint(value = a.value * b, bound = isnothing(a.bound) ? nothing : a.bound * b) +Base.:/(a::Constraint, b::Real) = + Constraint(value = a.value / b, bound = isnothing(a.bound) ? nothing : a.bound / b) """ $(TYPEDSIGNATURES) diff --git a/src/constraint_tree.jl b/src/constraint_tree.jl index c198430..f2c627a 100644 --- a/src/constraint_tree.jl +++ b/src/constraint_tree.jl @@ -186,10 +186,10 @@ single bound (of precise length 1) for creating all variables of the same constraint, or an iterable object of same length as `keys` with individual bounds for each variable in the same order as `keys`. -The individual bounds should be of type [`Bound`](@ref). To pass a single -interval bound for all variables, it is impossible to use a tuple (since its -length is 2); in such case use `bound = Ref((minimum, maximum))`, which has the -correct length. +The individual bounds should be subtypes of [`Bound`](@ref), or nothing. To +pass a single interval bound for all variables, it is impossible to use a tuple +(since its length is measured as 2); in such case use `bound = Ref((minimum, +maximum))`, which has the correct length. """ function variables(; keys::AbstractVector{Symbol}, bounds = nothing) bs = diff --git a/test/misc.jl b/test/misc.jl index c9dd568..0670ec3 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -33,9 +33,13 @@ end end @testset "Constraints" begin - @test C.bound(C.variable(bound = 123.0)) == 123.0 + @test C.bound(C.variable(bound = 123.0)).equal_to == 123.0 @test C.value(C.variable(bound = 123.0)).idxs == [1] - @test C.bound(2 * -convert(C.Constraint, (C.variable(bound = 123.0))) / 2) == -123.0 + @test C.bound(-convert(C.Constraint, (C.variable(bound = 123.0))) * 2 / 2).equal_to == + -123.0 + @test let x = C.bound(-convert(C.Constraint, (C.variable(bound = (-1, 2))) * 2 / 2)) + (x.lower, x.upper) == (1.0, -2.0) + end x = C.variable().value s = :a^C.Constraint(x, 5.0) + :b^C.Constraint(x * x - x, (4.0, 6.0)) From ac52726105f7c08da7989051e9cafb0ee212efbd Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Thu, 21 Dec 2023 11:45:54 +0100 Subject: [PATCH 2/2] try fixing the doc ref --- docs/src/metabolic-modeling.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/src/metabolic-modeling.jl b/docs/src/metabolic-modeling.jl index e1f9aa4..c76a11c 100644 --- a/docs/src/metabolic-modeling.jl +++ b/docs/src/metabolic-modeling.jl @@ -348,7 +348,8 @@ Dict(k => v.fluxes.R_BIOMASS_Ecoli_core_w_GAM for (k, v) in result.community) # You can thus, e.g., set a single bound: c.exchanges.oxygen.bound = C.Between(-20.0, 20.0) -# ...or rebuild a whole constraint (using a tuple shortcut for [`Between`](@ref)): +# ...or rebuild a whole constraint (using a tuple shortcut for +# [`ConstraintTrees.Between`](@ref)): c.exchanges.biomass = C.Constraint(c.exchanges.biomass.value, (-20, 20)) # ...or even add new constraints, here using the index syntax for demonstration: