Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

make the bounds type open to extension #15

Merged
merged 2 commits into from
Dec 21, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -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"
Expand Down
15 changes: 8 additions & 7 deletions docs/src/metabolic-modeling.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
exaexa marked this conversation as resolved.
Show resolved Hide resolved
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)
Expand Down Expand Up @@ -346,9 +346,10 @@ 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)
exaexa marked this conversation as resolved.
Show resolved Hide resolved

# ...or rebuild a whole constraint:
# ...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:
Expand Down
20 changes: 8 additions & 12 deletions docs/src/quadratic-optimization.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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.)
Expand All @@ -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),
),
)

Expand Down Expand Up @@ -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)
Expand Down
55 changes: 49 additions & 6 deletions src/bound.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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}
42 changes: 19 additions & 23 deletions src/constraint.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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`
Expand All @@ -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)
Expand Down
8 changes: 4 additions & 4 deletions src/constraint_tree.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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 =
Expand Down
8 changes: 6 additions & 2 deletions test/misc.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down