Merge pull request #4 from COBREXA/mk-qp

support affine and quadratic values&constraints

Project.toml

@@ -16,11 +16,12 @@ JuMP = "1"
 julia = "1.6"
 
 [extras]
+Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 SBML = "e5567a89-2604-4b09-9718-f5f78e97c3bb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Downloads", "GLPK", "JuMP", "SBML", "Test"]
+test = ["Clarabel", "Downloads", "GLPK", "JuMP", "SBML", "Test"]

docs/Project.toml

@@ -1,8 +1,12 @@
 [deps]
+Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 SBML = "e5567a89-2604-4b09-9718-f5f78e97c3bb"
+
+[compat]
+Documenter = "<1"

docs/src/index.md

@@ -3,5 +3,5 @@
 
 ```@autodocs
 Modules = [ConstraintTrees]
-Pages = ["ConstraintTrees.jl"]
+Pages = ["src/ConstraintTrees.jl"]
 ```

docs/src/metabolic-modeling.jl

@@ -25,7 +25,7 @@ ecoli = SBML.readSBML("e_coli_core.xml")
 # Let's first build the constrained representation of the problem. First, we
 # will need a variable for each of the reactions in the model.
 
-c = C.allocate_variables(keys = Symbol.(keys(ecoli.reactions)))
+c = C.variables(keys = Symbol.(keys(ecoli.reactions)))
 
 @test length(C.elems(c)) == length(ecoli.reactions) #src
 
@@ -71,7 +71,7 @@ c[:fluxes][:R_PFK]
 # variables to their valid bounds as defined by the model:
 rxn_constraints =
     let rxn_bounds = Symbol.(keys(ecoli.reactions)) .=> zip(SBML.flux_bounds(ecoli)...)
-        C.make_constraint_tree(
+        C.constraint_tree(
             r => C.Constraint(value = c.fluxes[r].value, bound = (lb, ub)) for
             (r, ((lb, _), (ub, _))) in rxn_bounds # SBML units are ignored for simplicity
         )
@@ -105,14 +105,33 @@ collect(keys(c))
 # values and making constraints.
 sum(C.value.(values(c.fluxes)))
 
+# ### Affine values
+#
+# To simplify various modeling goals (mainly calculation of various kinds of
+# "distances"), the values support inclusion of an affine element -- the
+# variable with index 0 is assumed to be the "affine unit", and its assigned
+# value is fixed at `1.0`.
+
+# To demonstrate, let's make a small system with 2 variables.
+system = C.variables(keys = [:x, :y])
+
+# To add an affine element to a `Value`, simply add it as a `Real`
+# number, as in the linear transformations below:
+system =
+    :original_coords^system *
+    :transformed_coords^C.constraint_tree(
+        :xt => C.Constraint(value = 1 + system.x.value + 4 + system.y.value),
+        :yt => C.Constraint(value = 0.1 * (3 - system.y.value)),
+    )
+
 # ## Adding combined constraints
 
 # Metabolic modeling relies on the fact that the total rates of any metabolite
 # getting created and consumed by the reaction equals to zero (which
 # corresponds to conservation of mass). We can now add corresponding
 # "stoichiometric" network constraints by following the reactants and products
 # in the SBML structure:
-stoi_constraints = C.make_constraint_tree(
+stoi_constraints = C.constraint_tree(
     Symbol(m) => C.Constraint(
         value = -sum(
             (
@@ -151,6 +170,35 @@ c *=
         ),
     );
 
+# ## Solution trees
+#
+# To aid exploration of variable assignments in the constraint trees, we can
+# convert them to *solution trees*. These have the very same structure as
+# constraint trees, but carry only the "solved" constraint values instead of
+# full constraints.
+#
+# Let's demonstrate this quickly on the example of `system` with affine
+# variables from above. First, let's assume that someone solved the system (in
+# some way) and produced a solution of variables as follows:
+solution = [1.0, 5.0] # corresponds to :x and :y in order.
+
+# Solution tree is constructed in a straightforward manner:
+st = C.solution_tree(system, solution)
+
+# We can now check the values of the original values
+(st.original_coords.x, st.original_coords.y)
+
+@test isapprox(st.original_coords.x, 1.0) #src
+@test isapprox(st.original_coords.y, 5.0) #src
+
+# The other constraints automatically get their values that correspond to the
+# overall variable assignment:
+st_ = st.transformed_coords;
+(st_.xt, st_.yt)
+
+@test isapprox(st_.xt, 11.0) #src
+@test isapprox(st_.yt, -0.2) #src
+
 # ## Solving the constraint system using JuMP
 #
 # We can make a small function that throws our model into JuMP, optimizes it,
@@ -225,17 +273,13 @@ c =
 
 # We can create additional variables that represent total community intake of
 # oxygen, and total community production of biomass:
-c +=
-    :exchanges^C.allocate_variables(
-        keys = [:oxygen, :biomass],
-        bounds = [(-10.0, 10.0), nothing],
-    )
+c += :exchanges^C.variables(keys = [:oxygen, :biomass], bounds = [(-10.0, 10.0), nothing])
 
 # These can be constrained so that the total influx (or outflux) of each of the
 # registered metabolites is in fact equal to total consumption or production by
 # each of the species:
 c *=
-    :exchange_constraints^C.make_constraint_tree(
+    :exchange_constraints^C.constraint_tree(
         :oxygen => C.Constraint(
             value = c.exchanges.oxygen.value - c.community.species1.fluxes.R_EX_o2_e.value -
                     c.community.species2.fluxes.R_EX_o2_e.value,

docs/src/quadratic-optimization.jl

@@ -0,0 +1,153 @@
+
+# # Example: Quadratic optimization
+#
+# In this example we demonstrate the use of quadratic constraints and values.
+# We assume that the reader is already familiar with the construction of
+# `ConstraintTree`s; if not, it is advisable to read the previous part
+# of the documentation first.
+#
+# In short, quadratic values and constraints are expressed similarly as other
+# contents of the constraint trees using types `QValue` and
+# `QConstraint`, which are quadratic alikes of the linear
+# `Value` and `Constraint`.
+#
+# ## Working with quadratic values and constraints
+#
+# Algebraically, you can construct `QValue`s simply by multiplying the linear
+# `Value`s:
+
+import ConstraintTrees as C
+
+system = C.variables(keys = [:x, :y, :z]);
+qv = system.x.value * (system.y.value + 2 * system.z.value)
+
+@test qv.idxs == [(1, 2), (1, 3)] #src
+@test qv.weights == [1.0, 2.0] #src
+
+# As with `Value`s, the `QValue`s can be easily combined, giving a nice way to
+# specify e.g. weighted sums of squared errors with respect to various
+# directions. Let's make a tiny helper first:
+squared(x) = x * x
+
+# Now, we can play with common representations of error values:
+error_val =
+    squared(system.x.value + system.y.value - 1) +
+    squared(system.y.value + 5 * system.z.value - 3)
+
+# 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.QConstraint(qvalue = error_val, bound = (0.0, 100.0))
+
+# Let's pretend someone has solved the system, and see how much "error" the
+# solution has:
+solution = [1.0, 2.0, -1.0];
+st = C.solution_tree(system, solution);
+st.error
+
+# ...not bad for a first guess.
+
+# ## Building quadratic systems
+#
+# Let's create a small quadratic system that finds the closest distance between
+# an ellipse and a line and let some of the conic solvers available in JuMP
+# solve it. First, let's make a representation of a point in 2D:
+point = C.variables(keys = [:x, :y])
+
+# We can create a small system that constraints the point to stay within a
+# simple elliptical area centered around `(0.0, 10.0)`:
+ellipse_system = C.constraint_tree(
+    :point => point,
+    :in_area => C.QConstraint(
+        qvalue = squared(point.x.value) / 4 + squared(10.0 - point.y.value),
+        bound = (-Inf, 1.0),
+    ),
+);
+
+# We now create another small system that constraints another point to stay on
+# a line that crosses `(0, 0)` and `(1, 1)`. We could do this using a
+# dot-product representation of line, but that would lead to issues later
+# (mainly, the solver that we are planning to use only supports positive
+# definite quadratic forms as constraints). Instead, let's use a
+# single-variable-parametrized line equation.
+line_param = C.variable().value;
+line_system =
+    :point^C.constraint_tree(
+        :x => C.Constraint(value = 0 + 1 * line_param),
+        :y => C.Constraint(value = 0 + 1 * line_param),
+    );
+
+# Finally, let's connect the systems using `+` operator and add the objective
+# that would minimize the distance of the points:
+s = :ellipse^ellipse_system + :line^line_system;
+
+s *=
+    :objective^C.QConstraint(
+        qvalue = squared(s.ellipse.point.x.value - s.line.point.x.value) +
+                 squared(s.ellipse.point.y.value - s.line.point.y.value),
+    );
+# (Note that if we used `*` to connect the systems, the variables from the
+# definition of `point` would not be duplicated, and various non-interesting
+# logic errors would follow.)
+
+# ## Solving quadratic systems with JuMP
+#
+# To solve the above system, we need a matching solver that can work with
+# quadratic constraints. Also, we need to create the function that translates
+# the constraints into JuMP `Model`s to support the quadratic constraints.
+import JuMP
+function optimized_vars(cs::C.ConstraintTree, objective::Union{C.Value,C.QValue}, optimizer)
+    model = JuMP.Model(optimizer)
+    JuMP.@variable(model, x[1:C.var_count(cs)])
+    if objective isa C.Value
+        JuMP.@objective(model, JuMP.MAX_SENSE, C.value_product(objective, x))
+    elseif objective isa C.QValue
+        JuMP.@objective(model, JuMP.MAX_SENSE, C.qvalue_product(objective, x))
+    end
+    function add_constraint(c::C.Constraint)
+        if c.bound isa Float64
+            JuMP.@constraint(model, C.value_product(c.value, x) == c.bound)
+        elseif c.bound isa Tuple{Float64,Float64}
+            val = C.value_product(c.value, x)
+            isinf(c.bound[1]) || JuMP.@constraint(model, val >= c.bound[1])
+            isinf(c.bound[2]) || JuMP.@constraint(model, val <= c.bound[2])
+        end
+    end
+    function add_constraint(c::C.QConstraint)
+        if c.bound isa Float64
+            JuMP.@constraint(model, C.qvalue_product(c.qvalue, x) == c.bound)
+        elseif c.bound isa Tuple{Float64,Float64}
+            val = C.qvalue_product(c.qvalue, x)
+            isinf(c.bound[1]) || JuMP.@constraint(model, val >= c.bound[1])
+            isinf(c.bound[2]) || JuMP.@constraint(model, val <= c.bound[2])
+        end
+    end
+    function add_constraint(c::C.ConstraintTree)
+        add_constraint.(values(c))
+    end
+    add_constraint(cs)
+    JuMP.set_silent(model)
+    JuMP.optimize!(model)
+    JuMP.value.(model[:x])
+end
+
+# We can now load a suitable optimizer and solve the system:
+import Clarabel
+st = C.solution_tree(s, optimized_vars(s, -s.objective.qvalue, Clarabel.Optimizer))
+
+# If the optimization worked well, we can nicely get out the position of the
+# closest point to the line that is in the elliptical area:
+(st.ellipse.point.x, st.ellipse.point.y)
+
+@test isapprox(st.ellipse.point.x, 1.7888553691812248, atol = 1e-3) #src
+@test isapprox(st.ellipse.point.y, 9.552787347840578, atol = 1e-3) #src
+
+# ...as well as the position on the line that is closest to the ellipse:
+C.elems(st.line.point)
+
+@test isapprox(st.line.point.x, st.line.point.y, atol = 1e-3) #src
+@test isapprox(st.line.point.x, 5.670821358510901, atol = 1e-3) #src
+
+# ...and, with a bit of extra math, the minimized distance:
+sqrt(st.objective)
+
+@test isapprox(sqrt(st.objective), 5.489928950781118, atol = 1e-3) #src

docs/src/reference.md

@@ -1,30 +1,53 @@
 
 # Reference
 
-## Values (linear combinations of variables)
+## Values
+
+### Linear and affine values
+
+```@autodocs
+Modules = [ConstraintTrees]
+Pages = ["src/value.jl"]
+```
+
+### Quadratic values
 
 ```@autodocs
 Modules = [ConstraintTrees]
-Pages = ["value.jl"]
+Pages = ["src/qvalue.jl"]
 ```
 
 ## Constraints
 
 ```@autodocs
 Modules = [ConstraintTrees]
-Pages = ["constraint.jl"]
+Pages = ["src/bound.jl"]
+```
+
+### Linear and affine constraints
+
+```@autodocs
+Modules = [ConstraintTrees]
+Pages = ["src/constraint.jl"]
+```
+
+### Quadratic constraints
+
+```@autodocs
+Modules = [ConstraintTrees]
+Pages = ["src/qconstraint.jl"]
 ```
 
 ## Constraint trees
 
 ```@autodocs
 Modules = [ConstraintTrees]
-Pages = ["constraint_tree.jl"]
+Pages = ["src/constraint_tree.jl"]
 ```
 
 ## Solution trees
 
 ```@autodocs
 Modules = [ConstraintTrees]
-Pages = ["solution.jl"]
+Pages = ["src/solution_tree.jl"]
 ```

src/ConstraintTrees.jl

@@ -6,8 +6,11 @@ module ConstraintTrees
 using DocStringExtensions
 
 include("value.jl")
+include("qvalue.jl")
+include("bound.jl")
 include("constraint.jl")
+include("qconstraint.jl")
 include("constraint_tree.jl")
-include("solution.jl")
+include("solution_tree.jl")
 
 end # module ConstraintTrees

src/bound.jl

@@ -0,0 +1,17 @@
+
+"""
+$(TYPEDEF)
+
+Convenience shortcut for "interval" bound; consisting of lower and upper bound
+value.
+"""
+const IntervalBound = Tuple{Float64,Float64}
+
+"""
+$(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.
+"""
+const Bound = Union{Nothing,Float64,IntervalBound}