january cleanup

Project.toml

@@ -1,7 +1,7 @@
 name = "ConstraintTrees"
 uuid = "5515826b-29c3-47a5-8849-8513ac836620"
 authors = ["The developers of ConstraintTrees.jl"]
-version = "0.8.1"
+version = "0.9"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"

docs/src/1-metabolic-modeling.jl

@@ -212,7 +212,7 @@ c *=
 solution = [1.0, 5.0] # corresponds to :x and :y in order given in `variables`
 
 # A value tree for this solution is constructed in a straightforward manner:
-st = C.constraint_values(system, solution)
+st = C.substitute_values(system, solution)
 
 # We can now check the values of the original coordinates
 st.original_coords
@@ -263,7 +263,7 @@ optimal_variable_assignment = optimized_vars(c, c.objective.value, GLPK.Optimize
 
 # To explore the solution more easily, we can make a tree with values that
 # correspond to ones in our constraint tree:
-result = C.constraint_values(c, optimal_variable_assignment)
+result = C.substitute_values(c, optimal_variable_assignment)
 
 result.fluxes.R_BIOMASS_Ecoli_core_w_GAM
 
@@ -273,11 +273,15 @@ result.fluxes.R_PFK
 
 # Sometimes it is unnecessary to recover the values for all constraints, so we
 # are better off selecting just the right subtree:
-C.constraint_values(c.fluxes, optimal_variable_assignment)
+C.substitute_values(c.fluxes, optimal_variable_assignment)
 
 #
 
-C.constraint_values(c.objective, optimal_variable_assignment)
+C.substitute_values(c.objective, optimal_variable_assignment)
+
+# We'll save the `result` for future use at the end of this example:
+
+result_single_organism = result
 
 # ## Combining and extending constraint systems
 #
@@ -322,7 +326,7 @@ c *=
 
 # Let's see how much biomass are the two species capable of producing together:
 result =
-    C.constraint_values(c, optimized_vars(c, c.exchanges.biomass.value, GLPK.Optimizer))
+    C.substitute_values(c, optimized_vars(c, c.exchanges.biomass.value, GLPK.Optimizer))
 result.exchanges
 
 # Finally, we can iterate over all species in the small community and see how
@@ -374,8 +378,52 @@ c[:exchanges][:production_is_zero] = C.Constraint(c.exchanges.biomass.value, 0)
 delete!(c.exchanges, :production_is_zero)
 
 # In the end, the flux optimization yields an expectably different result:
-result =
-    C.constraint_values(c, optimized_vars(c, c.exchanges.biomass.value, GLPK.Optimizer))
+result_with_more_oxygen =
+    C.substitute_values(c, optimized_vars(c, c.exchanges.biomass.value, GLPK.Optimizer))
 result.exchanges
 
-@test result.exchanges.oxygen < -19.0 #src
+@test result_with_more_oxygen.exchanges.oxygen < -19.0 #src
+
+# ## Combining trees
+#
+# ConstraintTrees.jl defines its own version of `zip` function that can apply a
+# function to the contents of several trees, "zipping" them over the same keys
+# in the structure. This is vaguely similar but otherwise not related to the
+# `zip` from Julia base (similarly, ConstraintTrees.jl have their own specific
+# `map`).
+#
+# In practice, this allows you to create combined trees with various nice
+# properties very quickly. For example, you can find how much the values have
+# changed between our two communities:
+
+C.zip((x, y) -> y - x, result, result_with_more_oxygen, Float64)
+
+# The result is again a `Tree`, with the contained type specified by the last
+# argument (`Float64` in this case). We can explore it right away as the other
+# result trees. Also, it is possible to call this kind of function using the
+# Julia `do` notation, making the syntax a bit neater:
+
+difference = C.zip(result, result_with_more_oxygen, Float64) do x, y
+    y - x
+end
+
+# Exploring the difference works as expected:
+
+difference.community.species1.fluxes
+
+# For convenience in special cases, `zip` is also overloaded for 3 arguments.
+# We can, for a completely artificial example, check if the absolute flux
+# change was bigger in the first or in the second organism in the community
+# when compared to the original single-organism flux (which we luckily saved
+# above):
+
+changes = C.zip(
+    result.community.species1,
+    result.community.species2,
+    result_single_organism,
+    Bool,
+) do s1, s2, orig
+    abs(s1 - orig) > abs(s2 - orig)
+end
+
+changes.fluxes

docs/src/2-quadratic-optimization.jl

@@ -54,7 +54,7 @@ system = :vars^system * :error^C.Constraint(error_val, C.Between(0, 100))
 # 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.constraint_values(system, solution)
+st = C.substitute_values(system, solution)
 st.error
 
 # ...not bad for a first guess.
@@ -135,7 +135,7 @@ end
 # We can now load a suitable optimizer and solve the system by maximizing the
 # negative squared error:
 import Clarabel
-st = C.constraint_values(s, quad_optimized_vars(s, -s.objective.value, Clarabel.Optimizer))
+st = C.substitute_values(s, quad_optimized_vars(s, -s.objective.value, 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:

docs/src/3-mixed-integer-optimization.jl

@@ -18,7 +18,7 @@
 # # Example: Mixed integer optimization (MILP)
 #
 # This example demonstrates the extension of `ConstraintTree` bounds structures
-# to accomodate new kinds of problems. In particular, we create a new kind of
+# to accommodate new kinds of problems. In particular, we create a new kind of
 # `Bound` that is restricting the value to be a full integer, and then solve a
 # geometric problem with that.
 
@@ -66,7 +66,7 @@ function jump_constraint(m, x, v::C.Value, b::C.Between)
 end
 
 # JuMP does not support direct integrality constraints, so we need to make a
-# small disgression with a slack variable:
+# small digression with a slack variable:
 function jump_constraint(m, x, v::C.Value, b::IntegerFromTo)
     var = JuMP.@variable(m, integer = true)
     JuMP.@constraint(m, var >= b.from)
@@ -101,7 +101,7 @@ dice_system *=
 
 # For solving, we use GLPK (it has MILP capabilities).
 import GLPK
-dices_thrown = C.constraint_values(
+dices_thrown = C.substitute_values(
     dice_system,
     milp_optimized_vars(
         dice_system,
@@ -121,7 +121,7 @@ dices_thrown = C.constraint_values(
 vars = C.variables(keys = [:a, :b, :c], bounds = IntegerFromTo(1, 100))
 
 # For simpliclty, we make a shortcut for "values" in all variables:
-v = C.tree_map(vars, C.value, C.Value)
+v = C.map(C.value, vars, C.Value)
 
 # With that shortcut, the constraint tree constructs quite easily:
 triangle_system =
@@ -134,7 +134,7 @@ triangle_system =
 # We will need a solver that supports both quadratic and integer optimization:
 import SCIP
 triangle_sides =
-    C.constraint_values(
+    C.substitute_values(
         triangle_system,
         milp_optimized_vars(
             triangle_system,

src/ConstraintTrees.jl

@@ -47,7 +47,7 @@ The package is structured as follows:
   -- this forms the basis of the "tidy" algebra of constraints.
 - A variable assignment, which is typically the "solution" for a given
   constraint tree, can be combined with a [`ConstraintTree`](@ref) to create a
-  "value tree" via [`constraint_values`](@ref), which enables browsing of the
+  "value tree" via [`substitute_values`](@ref), which enables browsing of the
   optimization results in the very same structure as the input
   [`ConstraintTree`](@ref).
 

src/constraint.jl

@@ -27,17 +27,17 @@ becomes easily accessible for inspection and building other constraints.
 # Fields
 $(TYPEDFIELDS)
 """
-Base.@kwdef mutable struct Constraint{V,B}
+Base.@kwdef mutable struct Constraint
     "A value (typically a [`LinearValue`](@ref) or a [`QuadraticValue`](@ref))
     that describes what the constraint constraints."
-    value::V
+    value::Value
     "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
+    bound::MaybeBound = nothing
 
-    function Constraint(v::T, b::U = nothing) where {T<:Value,U<:MaybeBound}
-        new{T,U}(v, b)
+    function Constraint(v::Value, b::MaybeBound = nothing)
+        new(v, b)
     end
 end
 

src/constraint_tree.jl

@@ -188,7 +188,8 @@ $(TYPEDSIGNATURES)
 Allocate a single unnamed variable, returning a Constraint with an optionally
 specified `bound`.
 """
-variable(; bound = nothing) = Constraint(value = LinearValue([1], [1.0]); bound)
+variable(; bound = nothing, idx = 1) =
+    Constraint(value = LinearValue(Int[idx], Float64[1.0]); bound)
 
 """
 $(TYPEDSIGNATURES)
@@ -211,11 +212,26 @@ function variables(; keys::AbstractVector{Symbol}, bounds = nothing)
         length(bounds) == length(keys) ? bounds :
         error("lengths of bounds and keys differ for allocated variables")
     ConstraintTree(
-        k => Constraint(value = LinearValue(Int[i], Float64[1.0]), bound = b) for
-        ((i, k), b) in zip(enumerate(keys), bs)
+        k => variable(idx = i, bound = b) for ((i, k), b) in Base.zip(enumerate(keys), bs)
     )
 end
 
+"""
+$(TYPEDSIGNATURES)
+
+Allocate a variable for each item in a constraint tree (or any other kind of
+tree) and return a [`ConstraintTree`](@ref) with variables bounded by the
+`makebound` function, which converts a given tree element's value into a bound
+for the corresponding variable.
+"""
+function variables_for(makebound, ts::Tree)
+    var_idx = 0
+    map(ts, Constraint) do x
+        var_idx += 1
+        variable(idx = var_idx, bound = makebound(x))
+    end
+end
+
 #
 # Transforming the constraint trees
 #
@@ -226,13 +242,18 @@ $(TYPEDSIGNATURES)
 Substitute variable values from `y` into the constraint tree's constraint's
 values, getting a tree of "solved" constraint values for the given variable
 assignment.
+
+The third argument forces the output type (it is forwarded to
+[`map`](@ref)). The type gets defaulted from `eltype(y)`.
 """
-constraint_values(x::ConstraintTree, y::Vector{Float64}) =
-    tree_map(x, c -> substitute(value(c), y), Float64)
+substitute_values(x::ConstraintTree, y::AbstractVector, ::Type{T} = eltype(y)) where {T} =
+    map(x, T) do c
+        substitute(value(c), y)
+    end
 
 """
 $(TYPEDSIGNATURES)
 
-Fallback for [`constraint_values`](@ref) for a single constraint.
+Fallback for [`substitute_values`](@ref) for a single constraint.
 """
-constraint_values(x::Constraint, y::Vector{Float64}) = substitute(value(x), y)
+substitute_values(x::Constraint, y::AbstractVector, _ = eltype(y)) = substitute(value(x), y)

src/linear_value.jl

@@ -21,9 +21,11 @@ $(TYPEDEF)
 A representation of a "value" in a linear constrained optimization problem. The
 value is an affine linear combination of several variables.
 
-`LinearValue`s can be combined additively and multiplied by real-number constants.
+`LinearValue`s can be combined additively and multiplied by real-number
+constants.
 
-Multiplying two `LinearValue`s yields a quadratic form (in a [`QuadraticValue`](@ref)).
+Multiplying two `LinearValue`s yields a quadratic form (in a
+[`QuadraticValue`](@ref)).
 
 # Fields
 $(TYPEDFIELDS)
@@ -59,42 +61,63 @@ Base.:*(a::Real, b::LinearValue) = b * a
 Base.:*(a::LinearValue, b::Real) = LinearValue(idxs = a.idxs, weights = b .* a.weights)
 Base.:/(a::LinearValue, b::Real) = LinearValue(idxs = a.idxs, weights = a.weights ./ b)
 
-function Base.:+(a::LinearValue, b::LinearValue)
+"""
+$(TYPEDSIGNATURES)
+
+Helper function for implementing [`LinearValue`](@ref)-like objects. Given
+"sparse" representations of linear combinations, it computes a "merged" linear
+combination of 2 values added together.
+
+Zeroes are not filtered out.
+"""
+function add_sparse_linear_combination(
+    a_idxs::Vector{Int},
+    a_weights::Vector{T},
+    b_idxs::Vector{Int},
+    b_weights::Vector{T},
+)::Tuple{Vector{Int},Vector{T}} where {T}
     r_idxs = Int[]
     r_weights = Float64[]
     ai = 1
-    ae = length(a.idxs)
+    ae = length(a_idxs)
     bi = 1
-    be = length(b.idxs)
+    be = length(b_idxs)
     while ai <= ae && bi <= be
-        if a.idxs[ai] < b.idxs[bi]
-            push!(r_idxs, a.idxs[ai])
-            push!(r_weights, a.weights[ai])
+        if a_idxs[ai] < b_idxs[bi]
+            push!(r_idxs, a_idxs[ai])
+            push!(r_weights, a_weights[ai])
             ai += 1
-        elseif a.idxs[ai] > b.idxs[bi]
-            push!(r_idxs, b.idxs[bi])
-            push!(r_weights, b.weights[bi])
+        elseif a_idxs[ai] > b_idxs[bi]
+            push!(r_idxs, b_idxs[bi])
+            push!(r_weights, b_weights[bi])
             bi += 1
-        else # a.idxs[ai] == b.idxs[bi] -- merge case
-            push!(r_idxs, a.idxs[ai])
-            push!(r_weights, a.weights[ai] + b.weights[bi])
+        else # a_idxs[ai] == b_idxs[bi] -- merge case
+            push!(r_idxs, a_idxs[ai])
+            push!(r_weights, a_weights[ai] + b_weights[bi])
             ai += 1
             bi += 1
         end
     end
     while ai <= ae
-        push!(r_idxs, a.idxs[ai])
-        push!(r_weights, a.weights[ai])
+        push!(r_idxs, a_idxs[ai])
+        push!(r_weights, a_weights[ai])
         ai += 1
     end
     while bi <= be
-        push!(r_idxs, b.idxs[bi])
-        push!(r_weights, b.weights[bi])
+        push!(r_idxs, b_idxs[bi])
+        push!(r_weights, b_weights[bi])
         bi += 1
     end
-    LinearValue(idxs = r_idxs, weights = r_weights)
+    return (r_idxs, r_weights)
 end
 
+Base.:+(a::LinearValue, b::LinearValue) =
+    let
+        (idxs, weights) =
+            add_sparse_linear_combination(a.idxs, a.weights, b.idxs, b.weights)
+        LinearValue(; idxs, weights)
+    end
+
 """
 $(TYPEDSIGNATURES)