Merge pull request #37 from COBREXA/mk-sizehints

add sizehints to lin+quad value adding, supply a pairwise reduce

.github/workflows/ci.yml

@@ -37,6 +37,7 @@ jobs:
       - uses: julia-actions/julia-runtest@latest
         continue-on-error: ${{ matrix.version == 'nightly' }}
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v3
+      - uses: codecov/codecov-action@v4
         with:
           file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}

Project.toml

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

docs/src/1-metabolic-modeling.jl

@@ -131,6 +131,15 @@ collect(keys(c))
 # values when making new constraints:
 sum(C.value.(values(c.fluxes)))
 
+# Notably, ConstraintTrees provide their own implementation of `sum` which
+# typically works faster when adding many `Value`s together. The basic interface and results are otherwise
+# the same as with the `sum` from Base:
+
+C.sum(C.value.(values(c.fluxes)))
+
+#md # !!! danger "`Base.sum` vs. `ConstraintTrees.sum`"
+#md #     Since the `sum` from Base package is usually implemented as a left fold, it does not behave optimally when the temporary sub-results grow during the computation (and thus their addition becomes gradually slower). In turn, using the `Base.sum` for summing up [`LinearValue`](@ref ConstraintTrees.LinearValue)s and [`QuadraticValue`](@ref ConstraintTrees.QuadraticValue)s may take time quadratic in the number of added items. [`sum`](@ref ConstraintTrees.sum) from ConstraintTrees uses a different addition order which reduces the amount of large items added together (implemented by ["pairwise" `preduce`](@ref ConstraintTrees.preduce)), and in works in almost-linear time in most cases.
+
 # ### Affine values
 #
 # To simplify various modeling goals (mainly calculation of various kinds of
@@ -159,13 +168,13 @@ system =
 # in the SBML structure:
 stoi_constraints = C.ConstraintTree(
     Symbol(m) => C.Constraint(
-        value = -sum(
+        value = -C.sum(
             (
                 sr.stoichiometry * c.fluxes[Symbol(rid)].value for
                 (rid, r) in ecoli.reactions for sr in r.reactants if sr.species == m
             ),
             init = zero(C.LinearValue), # sometimes the sums are empty
-        ) + sum(
+        ) + C.sum(
             (
                 sr.stoichiometry * c.fluxes[Symbol(rid)].value for
                 (rid, r) in ecoli.reactions for sr in r.products if sr.species == m
@@ -193,10 +202,11 @@ c = c * :stoichiometry^stoi_constraints
 # We can save that information into the constraint system immediately:
 c *=
     :objective^C.Constraint(
-        sum(
+        C.sum(
             c.fluxes[Symbol(rid)].value * coeff for
             (rid, coeff) in (keys(ecoli.reactions) .=> SBML.flux_objective(ecoli)) if
-            coeff != 0.0
+            coeff != 0.0;
+            init = 0.0,
         ),
     )
 
@@ -364,7 +374,7 @@ Dict(k => v.fluxes.R_BIOMASS_Ecoli_core_w_GAM for (k, v) in result.community)
 c.exchanges.oxygen.bound = C.Between(-20.0, 20.0)
 
 # ...or rebuild a whole constraint (using a tuple shortcut for
-# [`ConstraintTrees.Between`](@ref)):
+# [`Between`](@ref ConstraintTrees.Between)):
 c.exchanges.biomass = C.Constraint(c.exchanges.biomass.value, (-20, 20))
 
 # ...or even add new constraints, here using the index syntax for demonstration:

src/constraint_tree.jl

@@ -114,9 +114,7 @@ $(TYPEDSIGNATURES)
 Find the expected count of variables in a [`QuadraticValue`](@ref). (This is a
 O(1) operation, relying on the co-lexicographical ordering of indexes.)
 """
-var_count(x::QuadraticValue) = isempty(x.idxs) ? 0 : let (_, max) = last(x.idxs)
-    max
-end
+var_count(x::QuadraticValue) = isempty(x.idxs) ? 0 : last(last(x.idxs))
 
 """
 $(TYPEDSIGNATURES)

src/linear_value.jl

@@ -82,6 +82,10 @@ function add_sparse_linear_combination(
     ae = length(a_idxs)
     bi = 1
     be = length(b_idxs)
+
+    sizehint!(r_idxs, ae + be)
+    sizehint!(r_weights, ae + be)
+
     while ai <= ae && bi <= be
         if a_idxs[ai] < b_idxs[bi]
             push!(r_idxs, a_idxs[ai])

src/quadratic_value.jl

@@ -112,6 +112,9 @@ function add_sparse_quadratic_combination(
     bi = 1
     be = length(b_idxs)
 
+    sizehint!(r_idxs, ae + be)
+    sizehint!(r_weights, ae + be)
+
     while ai <= ae && bi <= be
         if colex_le(a_idxs[ai], b_idxs[bi])
             push!(r_idxs, a_idxs[ai])

src/tree.jl

@@ -155,7 +155,7 @@ function traverse(f, x)
         end
     go(x) = begin
         f(x)
-        nothing
+        return nothing
     end
 
     go(x)
@@ -174,7 +174,7 @@ function itraverse(f, x)
         end
     go(ix, x) = begin
         f(ix, x)
-        nothing
+        return nothing
     end
 
     go((), x)

src/value.jl

@@ -37,3 +37,60 @@ overload for the purpose of having [`substitute_values`](@ref) to run on both
 [`Constraint`](@ref)s and [`Value`](@ref)s.
 """
 substitute_values(x::Value, y::AbstractVector, _ = eltype(y)) = substitute(x, y)
+
+"""
+$(TYPEDSIGNATURES)
+
+An alternative of `Base.reduce` which does a "pairwise" reduction in the shape
+of a binary merge tree, like in mergesort. In general this is a little more
+complex, but if the reduced value "grows" with more elements added (such as
+when adding a lot of [`LinearValue`](@ref)s together), this is able to prevent
+a complexity explosion by postponing "large" reducing operations as much as
+possible.
+
+In the specific case with adding lots of [`LinearValue`](@ref)s and
+[`QuadraticValue`](@ref)s together, this effectively squashes the reduction
+complexity from something around `O(n^2)` to `O(n)` (with a little larger
+constant factor.
+"""
+function preduce(op, xs; init = zero(eltype(xs)), stack_type = eltype(xs))
+    # This works by simulating integer increment and carry to organize the
+    # additions in a (mildly begin-biased) tree. `used` stores the integer,
+    # `val` the associated values.
+    used = Vector{Bool}()
+    val = Vector{stack_type}()
+
+    for item in xs
+        idx = 1
+        while true
+            if idx > length(used)
+                push!(used, false)
+                push!(val, init)
+            end
+            used[idx] || break
+            item = op(item, val[idx]) # collect the bit and carry
+            used[idx] = false
+            idx += 1
+        end
+        val[idx] = item # hit a zero, no more carrying
+        used[idx] = true
+    end
+    # collect all used bits
+    item = init
+    for idx = 1:length(used)
+        if used[idx]
+            item = op(item, val[idx])
+        end
+    end
+    return item
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Alias for [`preduce`](@ref) that uses `+` as the operation.
+
+Not as versatile as the `sum` from Base, but much faster for growing values
+like [`LinearValue`](@ref)s and [`QuadraticValue`](@ref)s.
+"""
+sum(xs; init = zero(eltype(xs))) = preduce(+, xs; init)