Merge pull request #15 from COBREXA/mk-open-bounds

make the bounds type open to extension

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"

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,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)
 
-# ...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:

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)

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}

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)

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 =

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))