From 0cd9128a7b04bfa37556835b101592f0e1ed550b Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Thu, 16 May 2024 09:09:53 +0200 Subject: [PATCH 01/14] add sizehints to lin+quad value adding --- src/linear_value.jl | 4 ++++ src/quadratic_value.jl | 3 +++ 2 files changed, 7 insertions(+) diff --git a/src/linear_value.jl b/src/linear_value.jl index cf5f6c5..a6e3845 100644 --- a/src/linear_value.jl +++ b/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]) diff --git a/src/quadratic_value.jl b/src/quadratic_value.jl index 2ab1284..2a28bc6 100644 --- a/src/quadratic_value.jl +++ b/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]) From a051fc896932cc7b31853c52b72f04ebbcffa8f4 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Sun, 19 May 2024 12:55:45 +0200 Subject: [PATCH 02/14] implement a preliminary version of pairwise reduce --- src/value.jl | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/src/value.jl b/src/value.jl index e9776a1..f7439ab 100644 --- a/src/value.jl +++ b/src/value.jl @@ -37,3 +37,36 @@ 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) + # TODO improve type stability here (it's veeeery far from optimal). + # TODO find a way to smuggle this into mapreduce + up(::Nothing, _, i) = i + up(next::Tuple, l, i) = + let + (next1, i1) = down(next, l) + up(next1, l + 1, op(i, i1)) + end + down(::Nothing, _) = (nothing, init) + down(next::Tuple, l) = + l == 0 ? (iterate(xs, last(next)), first(next)) : + let (next1, x1) = down(next, l - 1), (next2, x2) = down(next1, l - 1) + (next2, op(x1, x2)) + end + up(iterate(xs), 0, init) +end From 7edb31fa13e058e57886be4fd41c0ca496e7d913 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 5 Jun 2024 20:33:34 +0200 Subject: [PATCH 03/14] use a stacky version of parallel reduce --- src/value.jl | 41 +++++++++++++++++++++++++++-------------- 1 file changed, 27 insertions(+), 14 deletions(-) diff --git a/src/value.jl b/src/value.jl index f7439ab..825196a 100644 --- a/src/value.jl +++ b/src/value.jl @@ -53,20 +53,33 @@ In the specific case with adding lots of [`LinearValue`](@ref)s and complexity from something around `O(n^2)` to `O(n)` (with a little larger constant factor. """ -function preduce(op, xs; init) - # TODO improve type stability here (it's veeeery far from optimal). - # TODO find a way to smuggle this into mapreduce - up(::Nothing, _, i) = i - up(next::Tuple, l, i) = - let - (next1, i1) = down(next, l) - up(next1, l + 1, op(i, i1)) +function preduce(op, xs; init = zero(eltype(xs))) + n = length(xs) + n == 0 && return init + + # 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. + stksize = sizeof(typeof(n)) * 8 - leading_zeros(n) + used = fill(false, stksize) + val = fill(init, stksize) + + for item in xs + idx = 1 + while used[idx] + item = op(item, val[idx]) # collect the bit and carry + used[idx] = false + idx += 1 end - down(::Nothing, _) = (nothing, init) - down(next::Tuple, l) = - l == 0 ? (iterate(xs, last(next)), first(next)) : - let (next1, x1) = down(next, l - 1), (next2, x2) = down(next1, l - 1) - (next2, op(x1, x2)) + val[idx] = item # hit a zero, no more carrying + used[idx] = true + end + # collect all used bits + item = init + for idx = 1:stksize + if used[idx] + item = op(item, val[idx]) end - up(iterate(xs), 0, init) + end + return item end From f05fec9aacec8cf68ac010b577231585ef9a0464 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 5 Jun 2024 20:41:38 +0200 Subject: [PATCH 04/14] add a useful alias --- src/value.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/value.jl b/src/value.jl index 825196a..2d06b47 100644 --- a/src/value.jl +++ b/src/value.jl @@ -83,3 +83,13 @@ function preduce(op, xs; init = zero(eltype(xs))) 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) From b0dc51fa45a35ee8dce0619e286ae3afc32a4a4b Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 5 Jun 2024 20:53:19 +0200 Subject: [PATCH 05/14] add docs and tests --- docs/src/1-metabolic-modeling.jl | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/docs/src/1-metabolic-modeling.jl b/docs/src/1-metabolic-modeling.jl index a31ff86..e9abdec 100644 --- a/docs/src/1-metabolic-modeling.jl +++ b/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,7 +202,7 @@ 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 @@ -364,7 +373,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: From 648d87617bab1bd9843253b2663ac11d3a8e50e2 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 5 Jun 2024 20:56:04 +0200 Subject: [PATCH 06/14] fmt --- src/value.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/value.jl b/src/value.jl index 2d06b47..57a08a5 100644 --- a/src/value.jl +++ b/src/value.jl @@ -92,4 +92,4 @@ 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) +sum(xs; init = zero(eltype(xs))) = preduce(+, xs; init) From 8e32fd42c44e7dc821e1a09cca6980b68e1ef6d1 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 5 Jun 2024 21:14:18 +0200 Subject: [PATCH 07/14] everyone should halt for a moment and enjoy the type stability while it lasts --- src/value.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/value.jl b/src/value.jl index 57a08a5..131f15a 100644 --- a/src/value.jl +++ b/src/value.jl @@ -53,7 +53,7 @@ In the specific case with adding lots of [`LinearValue`](@ref)s and complexity from something around `O(n^2)` to `O(n)` (with a little larger constant factor. """ -function preduce(op, xs; init = zero(eltype(xs))) +function preduce(op, xs; init = zero(eltype(xs)), stack_type = eltype(xs)) n = length(xs) n == 0 && return init @@ -62,7 +62,7 @@ function preduce(op, xs; init = zero(eltype(xs))) # `val` the associated values. stksize = sizeof(typeof(n)) * 8 - leading_zeros(n) used = fill(false, stksize) - val = fill(init, stksize) + val = stack_type[init for _ in used] for item in xs idx = 1 From abacf14c31d6c434d3f9a2ab83d4f7c50603c660 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 5 Jun 2024 21:39:59 +0200 Subject: [PATCH 08/14] ok now we can't even have length --- src/value.jl | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/value.jl b/src/value.jl index 131f15a..e5902bb 100644 --- a/src/value.jl +++ b/src/value.jl @@ -60,13 +60,17 @@ 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. - stksize = sizeof(typeof(n)) * 8 - leading_zeros(n) - used = fill(false, stksize) - val = stack_type[init for _ in used] + used = Vector{Bool}() + val = Vector{stack_type}() for item in xs idx = 1 - while used[idx] + 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 @@ -76,7 +80,7 @@ function preduce(op, xs; init = zero(eltype(xs)), stack_type = eltype(xs)) end # collect all used bits item = init - for idx = 1:stksize + for idx = 1:length(used) if used[idx] item = op(item, val[idx]) end From aaab275b0fc2c98c2756492d86a36353427da096 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 5 Jun 2024 21:42:10 +0200 Subject: [PATCH 09/14] efficiency confirmed, bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 557e8fe..dbd20d8 100644 --- a/Project.toml +++ b/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" From 99abd55d203e66596ae6745122c44de0a6cc7273 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 5 Jun 2024 21:53:54 +0200 Subject: [PATCH 10/14] oh wait I actually killed length --- src/value.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/value.jl b/src/value.jl index e5902bb..ff062ef 100644 --- a/src/value.jl +++ b/src/value.jl @@ -54,9 +54,6 @@ 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)) - n = length(xs) - n == 0 && return init - # 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. From e8e82ac2ad9f5809da8143db3bf4835b97e0d1fa Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 5 Jun 2024 22:30:26 +0200 Subject: [PATCH 11/14] generators. --- docs/src/1-metabolic-modeling.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/src/1-metabolic-modeling.jl b/docs/src/1-metabolic-modeling.jl index e9abdec..3820c3c 100644 --- a/docs/src/1-metabolic-modeling.jl +++ b/docs/src/1-metabolic-modeling.jl @@ -205,7 +205,8 @@ c *= 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, ), ) From e4839230224f0f218d1094ecc5d757478b5891c0 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 5 Jun 2024 22:42:25 +0200 Subject: [PATCH 12/14] fix ccov --- .github/workflows/ci.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 11ead99..97241e4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -40,3 +40,4 @@ jobs: - uses: codecov/codecov-action@v3 with: file: lcov.info + token: ${{ secrets.CODECOV_TOKEN }} From 6569540bd5cf2055995c41b735596a6f6a4454cd Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 5 Jun 2024 22:55:48 +0200 Subject: [PATCH 13/14] fix coverage on julia 1.10 --- src/constraint_tree.jl | 4 +--- src/tree.jl | 4 ++-- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/constraint_tree.jl b/src/constraint_tree.jl index cd00e87..1eb0360 100644 --- a/src/constraint_tree.jl +++ b/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) diff --git a/src/tree.jl b/src/tree.jl index 80d897a..fbdb8e4 100644 --- a/src/tree.jl +++ b/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) From 43a72980e11cf17e82a257d7003e30a5d4cf80d0 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 5 Jun 2024 22:59:24 +0200 Subject: [PATCH 14/14] upd ccov --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 97241e4..5b0f563 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -37,7 +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 }}