From 615a6f5ea908fa04534baedcb8af206dd6b37437 Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <miroslav.kratochvil@uni.lu>
Date: Tue, 17 Jan 2023 15:32:30 +0100
Subject: [PATCH 01/23] remove obsolete print_banner() from the SIF container

(now it prints an error, thanks @fjconejero for spotting)
---
 cobrexa.def | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cobrexa.def b/cobrexa.def
index 0ed492006..97815389d 100644
--- a/cobrexa.def
+++ b/cobrexa.def
@@ -14,7 +14,7 @@ From: cylon-x/lcsb-biocore/julia-base:latest
     cd project
     julia --project -e 'import Pkg; Pkg.add("Tulip"); Pkg.instantiate(); Pkg.build(); Pkg.precompile();'
 
-    echo 'using COBREXA; COBREXA._print_banner();' > /project/.startup.jl
+    echo 'using COBREXA;' > /project/.startup.jl
 
     chmod -R a+rX $JULIA_DEPOT_PATH
 

From 98d9218a39c27959f2ef9edf48b44316f52542d5 Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <miroslav.kratochvil@uni.lu>
Date: Thu, 19 Jan 2023 08:44:32 +0100
Subject: [PATCH 02/23] LICENSE file type was not markdown

(do not confuse the poor renderers)
---
 LICENSE.md => LICENSE | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename LICENSE.md => LICENSE (100%)

diff --git a/LICENSE.md b/LICENSE
similarity index 100%
rename from LICENSE.md
rename to LICENSE

From 0bec5f2189822d31b69e4e67997b8e55578ae402 Mon Sep 17 00:00:00 2001
From: CompatHelper Julia <compathelper_noreply@julialang.org>
Date: Wed, 25 Jan 2023 01:40:12 +0000
Subject: [PATCH 03/23] CompatHelper: bump compat for SBML to 1, (keep existing
 compat)

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 3d2d1f1d5..1a7e17485 100644
--- a/Project.toml
+++ b/Project.toml
@@ -32,7 +32,7 @@ JuMP = "1"
 MAT = "0.10"
 MacroTools = "0.5.6"
 OrderedCollections = "1.4"
-SBML = "~1.3"
+SBML = "~1.3, 1"
 StableRNGs = "1.0"
 Tulip = "0.7.0, 0.8.0, 0.9.2"
 julia = "1.5"

From 4a3051965ee6e74e07875af4abe988440ce5e619 Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <exa.exa@gmail.com>
Date: Wed, 25 Jan 2023 09:24:31 +0100
Subject: [PATCH 04/23] be more careful about SBML version

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 1a7e17485..792fa8c6a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -32,7 +32,7 @@ JuMP = "1"
 MAT = "0.10"
 MacroTools = "0.5.6"
 OrderedCollections = "1.4"
-SBML = "~1.3, 1"
+SBML = "~1.3, ~1.4"
 StableRNGs = "1.0"
 Tulip = "0.7.0, 0.8.0, 0.9.2"
 julia = "1.5"

From 6712526def47dcd368975a1bf506d79e438afe7d Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <miroslav.kratochvil@uni.lu>
Date: Thu, 26 Jan 2023 11:10:56 +0100
Subject: [PATCH 05/23] remove a config from local CI

---
 .gitlab-ci.yml | 17 -----------------
 1 file changed, 17 deletions(-)

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 2199e9255..4971d5432 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -88,11 +88,6 @@ variables:
     - Invoke-Expression $Env:ARTENOLIS_SOFT_PATH"\julia\"$Env:JULIA_VER"\bin\julia --inline=yes --check-bounds=yes --color=yes --project=@. -e 'import Pkg; Pkg.test(; coverage = true)'"
     - exit $LASTEXITCODE
 
-.global_env_win8: &global_env_win8
-  tags:
-    - windows8
-  <<: *global_env_win
-
 .global_env_win10: &global_env_win10
   tags:
     - windows10
@@ -153,12 +148,6 @@ linux:julia1.6:
 # Additional platform&environment compatibility tests
 #
 
-windows8:julia1.8:
-  stage: test-compat
-  <<: *global_trigger_compat_tests
-  <<: *global_julia18
-  <<: *global_env_win8
-
 windows10:julia1.8:
   stage: test-compat
   <<: *global_trigger_compat_tests
@@ -171,12 +160,6 @@ mac:julia1.8:
   <<: *global_julia18
   <<: *global_env_mac
 
-windows8:julia1.6:
-  stage: test-compat
-  <<: *global_trigger_compat_tests
-  <<: *global_julia16
-  <<: *global_env_win8
-
 windows10:julia1.6:
   stage: test-compat
   <<: *global_trigger_compat_tests

From f3f6dc42f8c3168d2766c625556714ed23783c14 Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <miroslav.kratochvil@uni.lu>
Date: Fri, 27 Jan 2023 14:54:54 +0100
Subject: [PATCH 06/23] add an actual GPA/GRR to DNF parser

---
 Project.toml                         |   2 +
 src/COBREXA.jl                       |   5 +-
 src/base/utils/gene_associations.jl  | 199 +++++++++++++++++----------
 test/base/utils/gene_associations.jl |  15 ++
 4 files changed, 144 insertions(+), 77 deletions(-)
 create mode 100644 test/base/utils/gene_associations.jl

diff --git a/Project.toml b/Project.toml
index 792fa8c6a..a64444102 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,6 +14,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+PikaParser = "3bbf5609-3e7b-44cd-8549-7c69f321e792"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SBML = "e5567a89-2604-4b09-9718-f5f78e97c3bb"
@@ -32,6 +33,7 @@ JuMP = "1"
 MAT = "0.10"
 MacroTools = "0.5.6"
 OrderedCollections = "1.4"
+PikaParser = "0.5"
 SBML = "~1.3, ~1.4"
 StableRNGs = "1.0"
 Tulip = "0.7.0, 0.8.0, 0.9.2"
diff --git a/src/COBREXA.jl b/src/COBREXA.jl
index 9e75888f9..b59d0b17a 100644
--- a/src/COBREXA.jl
+++ b/src/COBREXA.jl
@@ -27,21 +27,22 @@ module COBREXA
 
 using Distributed
 using DistributedData
+using DocStringExtensions
 using HDF5
 using JSON
 using JuMP
 using LinearAlgebra
-using MAT
 using MacroTools
+using MAT
 using OrderedCollections
 using Random
 using Serialization
 using SparseArrays
 using StableRNGs
 using Statistics
-using DocStringExtensions
 
 import Base: findfirst, getindex, show
+import PikaParser as PP
 import Pkg
 import SBML # conflict with Reaction struct name
 
diff --git a/src/base/utils/gene_associations.jl b/src/base/utils/gene_associations.jl
index 109317900..6e7d86441 100644
--- a/src/base/utils/gene_associations.jl
+++ b/src/base/utils/gene_associations.jl
@@ -2,21 +2,30 @@
 """
 $(TYPEDSIGNATURES)
 
-Parse `SBML.GeneProductAssociation` structure to the simpler GeneAssociation.
-The input must be (implicitly) in a positive DNF.
+Parse `SBML.GeneProductAssociation` structure and convert it to a strictly
+positive DNF [`GeneAssociation`](@ref). Negation (`SBML.GPANot`) is not
+supported.
 """
 function _parse_grr(gpa::SBML.GeneProductAssociation)::GeneAssociation
-    parse_ref(x) =
-        typeof(x) == SBML.GPARef ? [x.gene_product] :
-        begin
-            @_models_log @warn "Could not parse a part of gene association, ignoring: $x"
-            String[]
+
+    function fold_and(dnfs::Vector{Vector{Vector{String}}})::Vector{Vector{String}}
+        if isempty(dnfs)
+            [String[]]
+        else
+            unique(unique(String[l; r]) for l in dnfs[1] for r in fold_and(dnfs[2:end]))
         end
-    parse_and(x) =
-        typeof(x) == SBML.GPAAnd ? vcat([parse_and(i) for i in x.terms]...) : parse_ref(x)
-    parse_or(x) =
-        typeof(x) == SBML.GPAOr ? vcat([parse_or(i) for i in x.terms]...) : [parse_and(x)]
-    return parse_or(gpa)
+    end
+
+    dnf(x::SBML.GPARef) = [[x.gene_product]]
+    dnf(x::SBML.GPAOr) = collect(Set(vcat(dnf.(x.terms)...)))
+    dnf(x::SBML.GPAAnd) = fold_and(dnf.(x.terms))
+    dnf(x) = throw(
+        DomainError(
+            x,
+            "unsupported gene product association contents of type $(typeof(x))",
+        ),
+    )
+    return dnf(gpa)
 end
 
 """
@@ -49,69 +58,109 @@ julia> _parse_grr("(YIL010W and YLR043C) or (YIL010W and YGR209C)")
 _parse_grr(s::String)::Maybe{GeneAssociation} = _maybemap(_parse_grr, _parse_grr_to_sbml(s))
 
 """
-$(TYPEDSIGNATURES)
-
-Internal helper for parsing the string GRRs into SBML data structures. More
-general than [`_parse_grr`](@ref).
+PikaParser grammar for stringy GRR expressions.
 """
-function _parse_grr_to_sbml(str::String)::Maybe{SBML.GeneProductAssociation}
-    s = str
-    toks = String[]
-    m = Nothing
-    while !isnothing(
-        begin
-            m = match(r"( +|[a-zA-Z0-9_-]+|[^ a-zA-Z0-9_()-]+|[(]|[)])(.*)", s)
-        end,
-    )
-        tok = strip(m.captures[1])
-        !isempty(tok) && push!(toks, tok)
-        s = m.captures[2]
+const _grr_grammar = begin
+    # characters that typically form the identifiers
+    isident(x::Char) =
+        isletter(x) ||
+        isdigit(x) ||
+        x == '_' ||
+        x == '-' ||
+        x == ':' ||
+        x == '.' ||
+        x == '\'' ||
+        x == '[' ||
+        x == ']' ||
+        x == '\x03' # a very ugly exception for badly parsed MAT files
+
+    # scanner helpers
+    eat(p) = m -> begin
+        last = 0
+        for i in eachindex(m)
+            p(m[i]) || break
+            last = i
+        end
+        last
     end
 
-    fail() = throw(DomainError(str, "Could not parse GRR"))
-
-    # shunting yard
-    ops = Symbol[]
-    vals = SBML.GeneProductAssociation[]
-    fold(sym, op) =
-        while !isempty(ops) && last(ops) == sym
-            r = pop!(vals)
-            l = pop!(vals)
-            pop!(ops)
-            push!(vals, op([l, r]))
-        end
-    for tok in toks
-        if tok in ["and", "AND", "&", "&&"]
-            push!(ops, :and)
-        elseif tok in ["or", "OR", "|", "||"]
-            fold(:and, SBML.GPAAnd)
-            push!(ops, :or)
-        elseif tok == "("
-            push!(ops, :paren)
-        elseif tok == ")"
-            fold(:and, SBML.GPAAnd)
-            fold(:or, SBML.GPAOr)
-            if isempty(ops) || last(ops) != :paren
-                fail()
-            else
-                pop!(ops)
-            end
-        else
-            push!(vals, SBML.GPARef(tok))
-        end
+    # eat one of keywords
+    kws(w...) = m -> begin
+        last = eat(isident)(m)
+        m[begin:last] in w ? last : 0
     end
 
-    fold(:and, SBML.GPAAnd)
-    fold(:or, SBML.GPAOr)
+    PP.make_grammar(
+        [:expr],
+        PP.flatten(
+            Dict(
+                :space => PP.first(PP.scan(eat(isspace)), PP.epsilon),
+                :id => PP.scan(eat(isident)),
+                :orop =>
+                    PP.first(PP.tokens("||"), PP.token('|'), PP.scan(kws("OR", "or"))),
+                :andop => PP.first(
+                    PP.tokens("&&"),
+                    PP.token('&'),
+                    PP.scan(kws("AND", "and")),
+                ),
+                :expr => PP.seq(:space, :orexpr, :space, PP.end_of_input),
+                :orexpr => PP.first(
+                    :or => PP.seq(:andexpr, :space, :orop, :space, :orexpr),
+                    :andexpr,
+                ),
+                :andexpr => PP.first(
+                    :and => PP.seq(:baseexpr, :space, :andop, :space, :andexpr),
+                    :baseexpr,
+                ),
+                :baseexpr => PP.first(
+                    :id,
+                    :parenexpr => PP.seq(
+                        PP.token('('),
+                        :space,
+                        :orexpr,
+                        :space,
+                        PP.token(')'),
+                    ),
+                ),
+            ),
+            Char,
+        ),
+    )
+end
 
-    if !isempty(ops) || length(vals) > 1
-        fail()
-    end
+_grr_grammar_open(m, _) =
+    m.rule == :expr ? Bool[0, 1, 0, 0] :
+    m.rule == :parenexpr ? Bool[0, 0, 1, 0, 0] :
+    m.rule in [:or, :and] ? Bool[1, 0, 0, 0, 1] :
+    m.rule in [:andexpr, :orexpr, :notexpr, :baseexpr] ? Bool[1] :
+    (false for _ in m.submatches)
+
+_grr_grammar_fold(m, _, subvals) =
+    m.rule == :id ? SBML.GPARef(m.view) :
+    m.rule == :and ? SBML.GPAAnd([subvals[1], subvals[5]]) :
+    m.rule == :or ? SBML.GPAOr([subvals[1], subvals[5]]) :
+    m.rule == :parenexpr ? subvals[3] :
+    m.rule == :expr ? subvals[2] : isempty(subvals) ? nothing : subvals[1]
+
+"""
+$(TYPEDSIGNATURES)
 
-    if isempty(vals)
-        nothing
+Internal helper for parsing the string GRRs into SBML data structures. More
+general than [`_parse_grr`](@ref).
+"""
+function _parse_grr_to_sbml(str::String)::Maybe{SBML.GeneProductAssociation}
+    all(isspace, str) && return nothing
+    tree = PP.parse_lex(_grr_grammar, str)
+    match = PP.find_match_at!(tree, :expr, 1)
+    if match > 0
+        return PP.traverse_match(
+            tree,
+            match,
+            open = _grr_grammar_open,
+            fold = _grr_grammar_fold,
+        )
     else
-        first(vals)
+        throw(DomainError(str, "cannot parse GRR"))
     end
 end
 
@@ -124,14 +173,14 @@ string.
 # Example
 ```
 julia> _unparse_grr(String, [["YIL010W", "YLR043C"], ["YIL010W", "YGR209C"]])
-"(YIL010W and YLR043C) or (YIL010W and YGR209C)"
+"(YIL010W && YLR043C) || (YIL010W && YGR209C)"
 ```
 """
-function _unparse_grr(::Type{String}, grr::GeneAssociation)::String
-    grr_strings = String[]
-    for gr in grr
-        push!(grr_strings, "(" * join([g for g in gr], " and ") * ")")
-    end
-    grr_string = join(grr_strings, " or ")
-    return grr_string
+function _unparse_grr(
+    ::Type{String},
+    grr::GeneAssociation;
+    and = " && ",
+    or = " || ",
+)::String
+    return join(("(" * join(gr, and) * ")" for gr in grr), or)
 end
diff --git a/test/base/utils/gene_associations.jl b/test/base/utils/gene_associations.jl
new file mode 100644
index 000000000..469b885d0
--- /dev/null
+++ b/test/base/utils/gene_associations.jl
@@ -0,0 +1,15 @@
+
+@testset "GRR parsing" begin
+    @test sort(
+        sort.(
+            COBREXA._parse_grr("(αλφα OR βητα\x03) AND (R2-D2's_gene OR prefix:su[ff]ix)")
+        ),
+    ) == [
+        ["R2-D2's_gene", "αλφα"],
+        ["R2-D2's_gene", "βητα\x03"],
+        ["prefix:su[ff]ix", "αλφα"],
+        ["prefix:su[ff]ix", "βητα\x03"],
+    ]
+    @test_throws DomainError COBREXA._parse_grr("(")
+    @test isnothing(COBREXA._parse_grr(" "))
+end

From eed3038f657f26b7f2c9c141c2b4f2a0ba779d12 Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <miroslav.kratochvil@uni.lu>
Date: Fri, 27 Jan 2023 15:18:59 +0100
Subject: [PATCH 07/23] fix the tests to match the new state

---
 test/base/types/CoreModel.jl | 4 ++--
 test/base/types/Reaction.jl  | 2 +-
 test/io/io.jl                | 2 +-
 3 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/test/base/types/CoreModel.jl b/test/base/types/CoreModel.jl
index 61e9ba96a..7e782d172 100644
--- a/test/base/types/CoreModel.jl
+++ b/test/base/types/CoreModel.jl
@@ -14,6 +14,6 @@ end
     @test Set(reactions(cm)) == Set(reactions(sm))
     @test Set(reactions(cm)) == Set(reactions(cm2))
 
-    @test reaction_gene_association(sm, reactions(sm)[1]) ==
-          reaction_gene_association(cm, reactions(sm)[1])
+    @test sort(sort.(reaction_gene_association(sm, reactions(sm)[1]))) ==
+          sort(sort.(reaction_gene_association(cm, reactions(sm)[1])))
 end
diff --git a/test/base/types/Reaction.jl b/test/base/types/Reaction.jl
index 3fa5a1c07..fa1b3dd94 100644
--- a/test/base/types/Reaction.jl
+++ b/test/base/types/Reaction.jl
@@ -32,7 +32,7 @@
     @test all(
         contains.(
             sprint(show, MIME("text/plain"), r1),
-            ["r1", "100.0", "g1 and g2", "glycolysis", "blah", "biocyc"],
+            ["r1", "100.0", "g1 && g2", "glycolysis", "blah", "biocyc"],
         ),
     )
 
diff --git a/test/io/io.jl b/test/io/io.jl
index 4eb037f0b..6399f0583 100644
--- a/test/io/io.jl
+++ b/test/io/io.jl
@@ -22,5 +22,5 @@
     @test n_reactions(reconmodel) == 10600
     recon_grrs = [r.grr for (i, r) in reconmodel.reactions if !isnothing(r.grr)]
     @test length(recon_grrs) == 5938
-    @test sum(length.(recon_grrs)) == 13903
+    @test sum(length.(recon_grrs)) == 31504
 end

From 4c0306132e672d2dd6ffc3cb7bff28e7ba688a99 Mon Sep 17 00:00:00 2001
From: "St. Elmo Wilken" <stelmozors@gmail.com>
Date: Fri, 27 Jan 2023 16:44:40 +0100
Subject: [PATCH 08/23] fix ordering issue

---
 test/base/types/CoreModelCoupled.jl | 9 ++-------
 1 file changed, 2 insertions(+), 7 deletions(-)

diff --git a/test/base/types/CoreModelCoupled.jl b/test/base/types/CoreModelCoupled.jl
index d64ba1ae9..b16be2985 100644
--- a/test/base/types/CoreModelCoupled.jl
+++ b/test/base/types/CoreModelCoupled.jl
@@ -19,13 +19,8 @@ end
           reaction_gene_association(cm, reactions(sm)[1])
 
     @test reaction_gene_association_vec(cm)[1:3] == [
-        [["b3916"], ["b1723"]],
-        [
-            ["b0902", "b0903", "b2579"],
-            ["b0902", "b0903"],
-            ["b0902", "b3114"],
-            ["b3951", "b3952"],
-        ],
+        [["b1723"], ["b3916"]],
+        [["b0902", "b0903", "b2579"], ["b0902", "b0903"], ["b3951", "b3952"], ["b0902", "b3114"]],
         [["b4025"]],
     ]
 end

From cd589ce34574d565ec6f959fab921d8dc875ab34 Mon Sep 17 00:00:00 2001
From: "St. Elmo Wilken" <stelmozors@gmail.com>
Date: Fri, 27 Jan 2023 16:45:09 +0100
Subject: [PATCH 09/23] fix ordering issue in GRRs manually

---
 test/analysis/gecko.jl   | 18 ++++++++++++++++++
 test/analysis/smoment.jl | 18 ++++++++++++++++++
 2 files changed, 36 insertions(+)

diff --git a/test/analysis/gecko.jl b/test/analysis/gecko.jl
index eab3a7083..90c387f9f 100644
--- a/test/analysis/gecko.jl
+++ b/test/analysis/gecko.jl
@@ -1,6 +1,24 @@
 @testset "GECKO" begin
     model = load_model(StandardModel, model_paths["e_coli_core.json"])
 
+        # fix ordering
+        model.reactions["ATPS4r"].grr = [
+            ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735"],
+            ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735", "b3739"]
+        ]
+        # model.reactions["GLCpts"].grr = [
+        #     ["b2417", "b1101", "b2415", "b2416"],
+        #     ["b1817", "b1818", "b1819", "b2415", "b2416"],
+        #     ["b2417", "b1621", "b2415", "b2416"],
+        # ]
+    
+        model.reactions["GLCpts"].grr = [
+            ["b2417", "b1621", "b2415", "b2416"],
+            ["b1817", "b1818", "b1819", "b2415", "b2416"],
+            ["b2417", "b1101", "b2415", "b2416"],
+        ]
+    
+
     get_reaction_isozymes =
         rid ->
             haskey(ecoli_core_reaction_kcats, rid) ?
diff --git a/test/analysis/smoment.jl b/test/analysis/smoment.jl
index cba241217..decf9b43e 100644
--- a/test/analysis/smoment.jl
+++ b/test/analysis/smoment.jl
@@ -3,6 +3,24 @@
 
     get_gene_product_mass = gid -> get(ecoli_core_gene_product_masses, gid, 0.0)
 
+    # fix ordering
+    model.reactions["ATPS4r"].grr = [
+        ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735"],
+        ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735", "b3739"]
+    ]
+    # model.reactions["GLCpts"].grr = [
+    #     ["b2417", "b1101", "b2415", "b2416"],
+    #     ["b1817", "b1818", "b1819", "b2415", "b2416"],
+    #     ["b2417", "b1621", "b2415", "b2416"],
+    # ]
+
+    model.reactions["GLCpts"].grr = [
+        ["b2417", "b1621", "b2415", "b2416"],
+        ["b1817", "b1818", "b1819", "b2415", "b2416"],
+        ["b2417", "b1101", "b2415", "b2416"],
+    ]
+
+
     get_reaction_isozyme =
         rid ->
             haskey(ecoli_core_reaction_kcats, rid) ?

From 4581ec8e1c9a55fbdbef64751f15ce2cf289339e Mon Sep 17 00:00:00 2001
From: "St. Elmo Wilken" <stelmozors@gmail.com>
Date: Fri, 27 Jan 2023 17:03:01 +0100
Subject: [PATCH 10/23] change test values (order of input changed, hence this
 fix)

---
 test/analysis/gecko.jl   | 32 +++++++++++++-------------------
 test/analysis/smoment.jl |  8 +-------
 2 files changed, 14 insertions(+), 26 deletions(-)

diff --git a/test/analysis/gecko.jl b/test/analysis/gecko.jl
index 90c387f9f..a2fedbf02 100644
--- a/test/analysis/gecko.jl
+++ b/test/analysis/gecko.jl
@@ -1,30 +1,24 @@
 @testset "GECKO" begin
     model = load_model(StandardModel, model_paths["e_coli_core.json"])
 
-        # fix ordering
-        model.reactions["ATPS4r"].grr = [
-            ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735"],
-            ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735", "b3739"]
-        ]
-        # model.reactions["GLCpts"].grr = [
-        #     ["b2417", "b1101", "b2415", "b2416"],
-        #     ["b1817", "b1818", "b1819", "b2415", "b2416"],
-        #     ["b2417", "b1621", "b2415", "b2416"],
-        # ]
-    
-        model.reactions["GLCpts"].grr = [
-            ["b2417", "b1621", "b2415", "b2416"],
-            ["b1817", "b1818", "b1819", "b2415", "b2416"],
-            ["b2417", "b1101", "b2415", "b2416"],
-        ]
-    
+    # fix ordering
+    model.reactions["ATPS4r"].grr = [
+        ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735"],
+        ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735", "b3739"]
+    ]
+
+    model.reactions["GLCpts"].grr = [
+        ["b2417", "b1101", "b2415", "b2416"],
+        ["b1817", "b1818", "b1819", "b2415", "b2416"],
+        ["b2417", "b1621", "b2415", "b2416"],
+    ]
 
     get_reaction_isozymes =
         rid ->
             haskey(ecoli_core_reaction_kcats, rid) ?
             collect(
                 Isozyme(
-                    Dict(grr .=> ecoli_core_protein_stoichiometry[rid][i]),
+                    Dict(grr .=> ecoli_core_protein_stoichiometry[rid][i]), # need to reverse this
                     ecoli_core_reaction_kcats[rid][i]...,
                 ) for (i, grr) in enumerate(reaction_gene_association(model, rid))
             ) : Isozyme[]
@@ -59,7 +53,7 @@
 
     @test isapprox(
         rxn_fluxes["BIOMASS_Ecoli_core_w_GAM"],
-        0.812827846796761,
+        0.8128851171859212,
         atol = TEST_TOLERANCE,
     )
 
diff --git a/test/analysis/smoment.jl b/test/analysis/smoment.jl
index decf9b43e..62bc6154e 100644
--- a/test/analysis/smoment.jl
+++ b/test/analysis/smoment.jl
@@ -8,11 +8,6 @@
         ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735"],
         ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735", "b3739"]
     ]
-    # model.reactions["GLCpts"].grr = [
-    #     ["b2417", "b1101", "b2415", "b2416"],
-    #     ["b1817", "b1818", "b1819", "b2415", "b2416"],
-    #     ["b2417", "b1621", "b2415", "b2416"],
-    # ]
 
     model.reactions["GLCpts"].grr = [
         ["b2417", "b1621", "b2415", "b2416"],
@@ -20,7 +15,6 @@
         ["b2417", "b1101", "b2415", "b2416"],
     ]
 
-
     get_reaction_isozyme =
         rid ->
             haskey(ecoli_core_reaction_kcats, rid) ?
@@ -53,7 +47,7 @@
 
     @test isapprox(
         rxn_fluxes["BIOMASS_Ecoli_core_w_GAM"],
-        0.8907273630431708,
+        0.8907347602586123,
         atol = TEST_TOLERANCE,
     )
 end

From 7224f462c368f13b03eb66f9ea8a6d9cec9b73f4 Mon Sep 17 00:00:00 2001
From: exaexa <exaexa@users.noreply.github.com>
Date: Fri, 27 Jan 2023 16:25:18 +0000
Subject: [PATCH 11/23] automatic formatting

triggered by @exaexa on PR #732
---
 test/analysis/gecko.jl              | 2 +-
 test/analysis/smoment.jl            | 2 +-
 test/base/types/CoreModelCoupled.jl | 7 ++++++-
 3 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/test/analysis/gecko.jl b/test/analysis/gecko.jl
index a2fedbf02..a0d895295 100644
--- a/test/analysis/gecko.jl
+++ b/test/analysis/gecko.jl
@@ -4,7 +4,7 @@
     # fix ordering
     model.reactions["ATPS4r"].grr = [
         ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735"],
-        ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735", "b3739"]
+        ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735", "b3739"],
     ]
 
     model.reactions["GLCpts"].grr = [
diff --git a/test/analysis/smoment.jl b/test/analysis/smoment.jl
index 62bc6154e..9a8dd3f35 100644
--- a/test/analysis/smoment.jl
+++ b/test/analysis/smoment.jl
@@ -6,7 +6,7 @@
     # fix ordering
     model.reactions["ATPS4r"].grr = [
         ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735"],
-        ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735", "b3739"]
+        ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735", "b3739"],
     ]
 
     model.reactions["GLCpts"].grr = [
diff --git a/test/base/types/CoreModelCoupled.jl b/test/base/types/CoreModelCoupled.jl
index b16be2985..28c66775d 100644
--- a/test/base/types/CoreModelCoupled.jl
+++ b/test/base/types/CoreModelCoupled.jl
@@ -20,7 +20,12 @@ end
 
     @test reaction_gene_association_vec(cm)[1:3] == [
         [["b1723"], ["b3916"]],
-        [["b0902", "b0903", "b2579"], ["b0902", "b0903"], ["b3951", "b3952"], ["b0902", "b3114"]],
+        [
+            ["b0902", "b0903", "b2579"],
+            ["b0902", "b0903"],
+            ["b3951", "b3952"],
+            ["b0902", "b3114"],
+        ],
         [["b4025"]],
     ]
 end

From 478c03acdc2d288d23e695d2ffbe6fbc360238fe Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <exa.exa@gmail.com>
Date: Fri, 27 Jan 2023 20:06:15 +0100
Subject: [PATCH 12/23] remove some duct tape

---
 test/analysis/gecko.jl   | 16 +-------
 test/analysis/smoment.jl | 14 +------
 test/data_static.jl      | 84 +---------------------------------------
 3 files changed, 4 insertions(+), 110 deletions(-)

diff --git a/test/analysis/gecko.jl b/test/analysis/gecko.jl
index a0d895295..e38bdd3a8 100644
--- a/test/analysis/gecko.jl
+++ b/test/analysis/gecko.jl
@@ -1,24 +1,12 @@
 @testset "GECKO" begin
     model = load_model(StandardModel, model_paths["e_coli_core.json"])
 
-    # fix ordering
-    model.reactions["ATPS4r"].grr = [
-        ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735"],
-        ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735", "b3739"],
-    ]
-
-    model.reactions["GLCpts"].grr = [
-        ["b2417", "b1101", "b2415", "b2416"],
-        ["b1817", "b1818", "b1819", "b2415", "b2416"],
-        ["b2417", "b1621", "b2415", "b2416"],
-    ]
-
     get_reaction_isozymes =
         rid ->
             haskey(ecoli_core_reaction_kcats, rid) ?
             collect(
                 Isozyme(
-                    Dict(grr .=> ecoli_core_protein_stoichiometry[rid][i]), # need to reverse this
+                    Dict(grr .=> fill(1.0, size(grr))),
                     ecoli_core_reaction_kcats[rid][i]...,
                 ) for (i, grr) in enumerate(reaction_gene_association(model, rid))
             ) : Isozyme[]
@@ -53,7 +41,7 @@
 
     @test isapprox(
         rxn_fluxes["BIOMASS_Ecoli_core_w_GAM"],
-        0.8128851171859212,
+        0.8128427019072836,
         atol = TEST_TOLERANCE,
     )
 
diff --git a/test/analysis/smoment.jl b/test/analysis/smoment.jl
index 9a8dd3f35..d9bc9049f 100644
--- a/test/analysis/smoment.jl
+++ b/test/analysis/smoment.jl
@@ -3,25 +3,13 @@
 
     get_gene_product_mass = gid -> get(ecoli_core_gene_product_masses, gid, 0.0)
 
-    # fix ordering
-    model.reactions["ATPS4r"].grr = [
-        ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735"],
-        ["b3736", "b3737", "b3738", "b3731", "b3732", "b3733", "b3734", "b3735", "b3739"],
-    ]
-
-    model.reactions["GLCpts"].grr = [
-        ["b2417", "b1621", "b2415", "b2416"],
-        ["b1817", "b1818", "b1819", "b2415", "b2416"],
-        ["b2417", "b1101", "b2415", "b2416"],
-    ]
-
     get_reaction_isozyme =
         rid ->
             haskey(ecoli_core_reaction_kcats, rid) ?
             argmax(
                 smoment_isozyme_speed(get_gene_product_mass),
                 Isozyme(
-                    Dict(grr .=> ecoli_core_protein_stoichiometry[rid][i]),
+                    Dict(grr .=> fill(1.0, size(grr))),
                     ecoli_core_reaction_kcats[rid][i]...,
                 ) for (i, grr) in enumerate(reaction_gene_association(model, rid))
             ) : nothing
diff --git a/test/data_static.jl b/test/data_static.jl
index a80d46d16..7af7a75ad 100644
--- a/test/data_static.jl
+++ b/test/data_static.jl
@@ -300,93 +300,11 @@ const ecoli_core_gene_product_masses = Dict{String,Float64}(
     "b2279" => 10.845,
 )
 
-const ecoli_core_protein_stoichiometry = Dict{String,Vector{Vector{Float64}}}(
-    #=
-    Data made up, each isozyme is assumed to be composed of
-    only one subunit each.
-    =#
-    "ACALD" => [[1.0], [1.0]],
-    "PTAr" => [[1.0], [1.0]],
-    "ALCD2x" => [[1.0], [1.0], [1.0]],
-    "PDH" => [[1.0, 1.0, 1.0]],
-    "PYK" => [[1.0], [1.0]],
-    "CO2t" => [[1.0]],
-    "MALt2_2" => [[1.0]],
-    "CS" => [[1.0]],
-    "PGM" => [[1.0], [1.0], [1.0]],
-    "TKT1" => [[1.0], [1.0]],
-    "ACONTa" => [[1.0], [1.0]],
-    "GLNS" => [[1.0], [1.0]],
-    "ICL" => [[1.0]],
-    "FBA" => [[1.0], [1.0], [1.0]],
-    "FORt2" => [[1.0], [1.0]],
-    "G6PDH2r" => [[1.0]],
-    "AKGDH" => [[1.0, 1.0, 1.0]],
-    "TKT2" => [[1.0], [1.0]],
-    "FRD7" => [[1.0, 1.0, 1.0, 1.0]],
-    "SUCOAS" => [[1.0, 1.0]],
-    "FBP" => [[1.0], [1.0]],
-    "ICDHyr" => [[1.0]],
-    "AKGt2r" => [[1.0]],
-    "GLUSy" => [[1.0, 1.0]],
-    "TPI" => [[1.0]],
-    "FORt" => [[1.0], [1.0]],
-    "ACONTb" => [[1.0], [1.0]],
-    "GLNabc" => [[1.0, 1.0, 1.0]],
-    "RPE" => [[1.0], [1.0]],
-    "ACKr" => [[1.0], [1.0], [1.0]],
-    "THD2" => [[1.0, 1.0]],
-    "PFL" => [[1.0, 1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0]],
-    "RPI" => [[1.0], [1.0]],
-    "D_LACt2" => [[1.0], [1.0]],
-    "TALA" => [[1.0], [1.0]],
-    "PPCK" => [[1.0]],
-    "ACt2r" => [[1.0]],
-    "NH4t" => [[1.0], [1.0]],
-    "PGL" => [[1.0]],
-    "NADTRHD" => [[1.0], [1.0, 1.0]],
-    "PGK" => [[1.0]],
-    "LDH_D" => [[1.0], [1.0]],
-    "ME1" => [[1.0]],
-    "PIt2r" => [[1.0], [1.0]],
-    "ATPS4r" => [
-        [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
-        [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
-    ],
-    "GLCpts" => [[1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0]],
-    "GLUDy" => [[1.0]],
-    "CYTBD" => [[1.0, 1.0], [1.0, 1.0]],
-    "FUMt2_2" => [[1.0]],
-    "FRUpts2" => [[1.0, 1.0, 1.0, 1.0, 1.0]],
-    "GAPD" => [[1.0]],
-    "H2Ot" => [[1.0], [1.0]],
-    "PPC" => [[1.0]],
-    "NADH16" => [[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]],
-    "PFK" => [[1.0], [1.0]],
-    "MDH" => [[1.0]],
-    "PGI" => [[1.0]],
-    "O2t" => [[1.0]],
-    "ME2" => [[1.0]],
-    "GND" => [[1.0]],
-    "SUCCt2_2" => [[1.0]],
-    "GLUN" => [[1.0], [1.0], [1.0]],
-    "ETOHt2r" => [[1.0]],
-    "ADK1" => [[1.0]],
-    "ACALDt" => [[1.0]],
-    "SUCDi" => [[1.0, 1.0, 1.0, 1.0]],
-    "ENO" => [[1.0]],
-    "MALS" => [[1.0], [1.0]],
-    "GLUt2r" => [[1.0]],
-    "PPS" => [[1.0]],
-    "FUM" => [[1.0], [1.0], [1.0]],
-)
-
 const ecoli_core_reaction_kcats = Dict{String,Vector{Tuple{Float64,Float64}}}(
     #=
     Data taken from Heckmann, David, et al. "Machine learning applied to enzyme
     turnover numbers reveals protein structural correlates and improves metabolic
-    models." Nature communications 9.1 (2018): 1-10. Assume forward and reverse
-    kcats are the same, and each isozyme has the same kcat.
+    models." Nature communications 9.1 (2018): 1-10.
     =#
     "ACALD" =>
         [(568.1130792316333, 568.1130792316333), (568.856126503717, 568.856126503717)],

From 74830b2d492e445e3512f85281c1b894135667af Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <exa.exa@gmail.com>
Date: Sat, 28 Jan 2023 15:33:24 +0100
Subject: [PATCH 13/23] make the GPA parsing completely independent of
 hashtable order

---
 src/base/utils/gene_associations.jl | 28 +++++++++++++++++++++++++++-
 test/analysis/gecko.jl              |  2 +-
 test/analysis/smoment.jl            |  2 +-
 test/io/io.jl                       |  2 +-
 4 files changed, 30 insertions(+), 4 deletions(-)

diff --git a/src/base/utils/gene_associations.jl b/src/base/utils/gene_associations.jl
index 6e7d86441..280fdb227 100644
--- a/src/base/utils/gene_associations.jl
+++ b/src/base/utils/gene_associations.jl
@@ -2,6 +2,30 @@
 """
 $(TYPEDSIGNATURES)
 
+A helper for producing predictable unique sequences. Might be faster if
+compacting would be done directly in sort().
+"""
+function _sortunique(x)
+    o = collect(x)
+    sort!(o)
+    put = prevind(o, firstindex(o))
+    for i in eachindex(o)
+        if put >= firstindex(o) && o[i] == o[put]
+            # we already have this one
+            continue
+        else
+            put = nextind(o, put)
+            if put != i
+                o[put] = o[i]
+            end
+        end
+    end
+    o[begin:put]
+end
+
+"""
+$(TYPEDSIGNATURES)
+
 Parse `SBML.GeneProductAssociation` structure and convert it to a strictly
 positive DNF [`GeneAssociation`](@ref). Negation (`SBML.GPANot`) is not
 supported.
@@ -12,7 +36,9 @@ function _parse_grr(gpa::SBML.GeneProductAssociation)::GeneAssociation
         if isempty(dnfs)
             [String[]]
         else
-            unique(unique(String[l; r]) for l in dnfs[1] for r in fold_and(dnfs[2:end]))
+            _sortunique(
+                _sortunique(String[l; r]) for l in dnfs[1] for r in fold_and(dnfs[2:end])
+            )
         end
     end
 
diff --git a/test/analysis/gecko.jl b/test/analysis/gecko.jl
index e38bdd3a8..62adb342b 100644
--- a/test/analysis/gecko.jl
+++ b/test/analysis/gecko.jl
@@ -41,7 +41,7 @@
 
     @test isapprox(
         rxn_fluxes["BIOMASS_Ecoli_core_w_GAM"],
-        0.8128427019072836,
+        0.8128851171859416,
         atol = TEST_TOLERANCE,
     )
 
diff --git a/test/analysis/smoment.jl b/test/analysis/smoment.jl
index d9bc9049f..a7b8a8068 100644
--- a/test/analysis/smoment.jl
+++ b/test/analysis/smoment.jl
@@ -35,7 +35,7 @@
 
     @test isapprox(
         rxn_fluxes["BIOMASS_Ecoli_core_w_GAM"],
-        0.8907347602586123,
+        0.8907680040632139,
         atol = TEST_TOLERANCE,
     )
 end
diff --git a/test/io/io.jl b/test/io/io.jl
index 6399f0583..f0e554906 100644
--- a/test/io/io.jl
+++ b/test/io/io.jl
@@ -22,5 +22,5 @@
     @test n_reactions(reconmodel) == 10600
     recon_grrs = [r.grr for (i, r) in reconmodel.reactions if !isnothing(r.grr)]
     @test length(recon_grrs) == 5938
-    @test sum(length.(recon_grrs)) == 31504
+    @test sum(length.(recon_grrs)) == 31503
 end

From b2a81cc59b7a403224336dbbf74bfb7ee08620ab Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <exa.exa@gmail.com>
Date: Sat, 28 Jan 2023 20:03:44 +0100
Subject: [PATCH 14/23] small fix 1

---
 src/base/types/StandardModel.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/base/types/StandardModel.jl b/src/base/types/StandardModel.jl
index 211a1e088..1d2cc88ae 100644
--- a/src/base/types/StandardModel.jl
+++ b/src/base/types/StandardModel.jl
@@ -187,7 +187,7 @@ Return the gene reaction rule in string format for reaction with `id` in `model`
 Return `nothing` if not available.
 """
 reaction_gene_association(model::StandardModel, id::String)::Maybe{GeneAssociation} =
-    _maybemap(identity, model.reactions[id].grr)
+    model.reactions[id].grr
 
 """
 $(TYPEDSIGNATURES)

From c47427c91984c69e590b523b54d18b6a38adbe7e Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <exa.exa@gmail.com>
Date: Sat, 28 Jan 2023 20:04:12 +0100
Subject: [PATCH 15/23] fix gra order in tests

---
 test/base/types/CoreModelCoupled.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/base/types/CoreModelCoupled.jl b/test/base/types/CoreModelCoupled.jl
index 28c66775d..2c2647a5a 100644
--- a/test/base/types/CoreModelCoupled.jl
+++ b/test/base/types/CoreModelCoupled.jl
@@ -21,10 +21,10 @@ end
     @test reaction_gene_association_vec(cm)[1:3] == [
         [["b1723"], ["b3916"]],
         [
-            ["b0902", "b0903", "b2579"],
             ["b0902", "b0903"],
-            ["b3951", "b3952"],
+            ["b0902", "b0903", "b2579"],
             ["b0902", "b3114"],
+            ["b3951", "b3952"],
         ],
         [["b4025"]],
     ]

From dcc348d74a0de8c042b975653b8187cd1662c511 Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <exa.exa@gmail.com>
Date: Sat, 28 Jan 2023 20:32:35 +0100
Subject: [PATCH 16/23] never going to use Set for compaction again

> Originally, Seth was a sky god, lord of the desert, master of storms,
> disorder, and warfare -- in general, a trickster. Set embodied the necessary
> and creative element of violence and disorder within the ordered world.
---
 src/base/utils/gene_associations.jl | 2 +-
 test/analysis/gecko.jl              | 2 +-
 test/analysis/smoment.jl            | 2 +-
 test/base/types/CoreModelCoupled.jl | 6 ++++--
 4 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/src/base/utils/gene_associations.jl b/src/base/utils/gene_associations.jl
index 280fdb227..de4532eaf 100644
--- a/src/base/utils/gene_associations.jl
+++ b/src/base/utils/gene_associations.jl
@@ -43,7 +43,7 @@ function _parse_grr(gpa::SBML.GeneProductAssociation)::GeneAssociation
     end
 
     dnf(x::SBML.GPARef) = [[x.gene_product]]
-    dnf(x::SBML.GPAOr) = collect(Set(vcat(dnf.(x.terms)...)))
+    dnf(x::SBML.GPAOr) = _sortunique(vcat(dnf.(x.terms)...))
     dnf(x::SBML.GPAAnd) = fold_and(dnf.(x.terms))
     dnf(x) = throw(
         DomainError(
diff --git a/test/analysis/gecko.jl b/test/analysis/gecko.jl
index 62adb342b..641b099bf 100644
--- a/test/analysis/gecko.jl
+++ b/test/analysis/gecko.jl
@@ -41,7 +41,7 @@
 
     @test isapprox(
         rxn_fluxes["BIOMASS_Ecoli_core_w_GAM"],
-        0.8128851171859416,
+        0.8129179015245396,
         atol = TEST_TOLERANCE,
     )
 
diff --git a/test/analysis/smoment.jl b/test/analysis/smoment.jl
index a7b8a8068..122f95a65 100644
--- a/test/analysis/smoment.jl
+++ b/test/analysis/smoment.jl
@@ -35,7 +35,7 @@
 
     @test isapprox(
         rxn_fluxes["BIOMASS_Ecoli_core_w_GAM"],
-        0.8907680040632139,
+        0.8907796509479097,
         atol = TEST_TOLERANCE,
     )
 end
diff --git a/test/base/types/CoreModelCoupled.jl b/test/base/types/CoreModelCoupled.jl
index 2c2647a5a..3bce19b95 100644
--- a/test/base/types/CoreModelCoupled.jl
+++ b/test/base/types/CoreModelCoupled.jl
@@ -15,8 +15,10 @@ end
     @test Set(reactions(cm)) == Set(reactions(sm))
     @test Set(reactions(cm)) == Set(reactions(cm2))
 
-    @test reaction_gene_association(sm, reactions(sm)[1]) ==
-          reaction_gene_association(cm, reactions(sm)[1])
+    @test all(
+        reaction_gene_association.(Ref(sm), reactions(sm)) .==
+        reaction_gene_association.(Ref(cm), reactions(sm)),
+    )
 
     @test reaction_gene_association_vec(cm)[1:3] == [
         [["b1723"], ["b3916"]],

From aeaf923052350e01413a80ad38cec353ebfecfba Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <miroslav.kratochvil@uni.lu>
Date: Fri, 10 Feb 2023 09:53:06 +0100
Subject: [PATCH 17/23] version bump for 1.4.4

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index a64444102..9f5b2465a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "COBREXA"
 uuid = "babc4406-5200-4a30-9033-bf5ae714c842"
 authors = ["The developers of COBREXA.jl"]
-version = "1.4.3"
+version = "1.4.4"
 
 [deps]
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"

From a13f98d85694ea462e4abe50dca30861fd827a05 Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <exa.exa@gmail.com>
Date: Sun, 12 Feb 2023 12:06:21 +0100
Subject: [PATCH 18/23] clean up SBML unit-error throwing

---
 src/base/types/SBMLModel.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/base/types/SBMLModel.jl b/src/base/types/SBMLModel.jl
index e07e7a31b..96c0633b8 100644
--- a/src/base/types/SBMLModel.jl
+++ b/src/base/types/SBMLModel.jl
@@ -139,7 +139,7 @@ function bounds(model::SBMLModel)::Tuple{Vector{Float64},Vector{Float64}}
                 if unit != common_unit
                     throw(
                         DomainError(
-                            units_in_sbml,
+                            unit,
                             "The SBML file uses multiple units; loading would need conversion",
                         ),
                     )

From 64023b3a2c6268a210528e58a5a88b5e117926e0 Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <miroslav.kratochvil@uni.lu>
Date: Wed, 3 May 2023 10:43:38 +0200
Subject: [PATCH 19/23] implement expression-limited fluxes (E-Flux like
 algorihm)

Closes #552
---
 src/analysis/expressions.jl                   | 11 +++++
 .../types/wrappers/ExpressionLimitedModel.jl  | 48 +++++++++++++++++++
 src/base/utils/eflux.jl                       | 16 +++++++
 src/reconstruction/expresisons.jl             | 10 ++++
 test/analysis/expressions.jl                  | 16 +++++++
 5 files changed, 101 insertions(+)
 create mode 100644 src/analysis/expressions.jl
 create mode 100644 src/base/types/wrappers/ExpressionLimitedModel.jl
 create mode 100644 src/base/utils/eflux.jl
 create mode 100644 src/reconstruction/expresisons.jl
 create mode 100644 test/analysis/expressions.jl

diff --git a/src/analysis/expressions.jl b/src/analysis/expressions.jl
new file mode 100644
index 000000000..32e82d37e
--- /dev/null
+++ b/src/analysis/expressions.jl
@@ -0,0 +1,11 @@
+
+"""
+$(TYPEDSIGNATURES)
+
+Build an [`ExpressionLimitedModel`](@ref).
+"""
+make_expression_limited_model(
+    model::MetabolicModel;
+    relative_expression::Dict{String,Float64},
+    bounding_function::Function = expression_probabilistic_bounds,
+) = ExpressionLimitedModel(; relative_expression, bounding_function, inner = model)
diff --git a/src/base/types/wrappers/ExpressionLimitedModel.jl b/src/base/types/wrappers/ExpressionLimitedModel.jl
new file mode 100644
index 000000000..e4f30612e
--- /dev/null
+++ b/src/base/types/wrappers/ExpressionLimitedModel.jl
@@ -0,0 +1,48 @@
+
+"""
+    $(TYPEDEF)
+
+[`ExpressionLimitedModel`](@ref) follows the methodology of the E-Flux algorithm to
+constraint the flux through the reactions in order to simulate the limited
+expression of genes (and thus the limited availability of the gene products).
+
+Use [`make_expression_limited_model`](@ref) or [`with_expression_limits`](@ref)
+to construct the models.
+
+E-Flux algorithm is closer described by: *Colijn, Caroline, Aaron Brandes,
+Jeremy Zucker, Desmond S. Lun, Brian Weiner, Maha R. Farhat, Tan-Yun Cheng, D.
+Branch Moody, Megan Murray, and James E. Galagan. "Interpreting expression data
+with metabolic flux models: predicting Mycobacterium tuberculosis mycolic acid
+production." PLoS computational biology 5, no. 8 (2009): e1000489*.
+
+# Fields
+$(TYPEDFIELDS)
+"""
+Base.@kwdef struct ExpressionLimitedModel <: ModelWrapper
+    """
+    Relative gene expression w.r.t. to some chosen reference; the normalization
+    and scale of the values should match the expectations of the
+    `bounding_function`.
+    """
+    relative_expression::Dict{String,Float64}
+
+    "The wrapped model."
+    inner::MetabolicModel
+
+    """
+    Function used to calculate the new reaction bounds from expression.
+    In [`make_expression_limited_model`](@ref) this defaults to
+    [`expression_probabilistic_bounds`](@ref).
+    """
+    bounding_function::Function
+end
+
+COBREXA.unwrap_model(m::ExpressionLimitedModel) = m.inner
+
+function COBREXA.bounds(m::ExpressionLimitedModel)::Tuple{Vector{Float64},Vector{Float64}}
+    (lbs, ubs) = bounds(m.inner)
+    lims = collect(
+        m.bounding_function(m.relative_expression, reaction_gene_association(m.inner, rid)) for rid in reactions(m.inner)
+    )
+    (lbs .* lims, ubs .* lims)
+end
diff --git a/src/base/utils/eflux.jl b/src/base/utils/eflux.jl
new file mode 100644
index 000000000..f5950c6fd
--- /dev/null
+++ b/src/base/utils/eflux.jl
@@ -0,0 +1,16 @@
+
+"""
+$(TYPEDSIGNATURES)
+
+Create E-Flux-like bounds for a reaction that requires gene products given by
+`grr`, with expression "choke" data given by relative probabilities (between 0
+and 1) in `relative_expression`.
+"""
+function expression_probabilistic_bounds(
+    relative_expression::Dict{String,Float64},
+    grr::Maybe{Vector{Vector{String}}},
+)::Float64
+    isnothing(grr) && return 1.0
+    lup(g) = get(relative_expression, g, 1.0)
+    1 - prod(1 - prod(lup.(gr)) for gr in grr)
+end
diff --git a/src/reconstruction/expresisons.jl b/src/reconstruction/expresisons.jl
new file mode 100644
index 000000000..c2bb80d2b
--- /dev/null
+++ b/src/reconstruction/expresisons.jl
@@ -0,0 +1,10 @@
+"""
+$(TYPEDSIGNATURES)
+
+Specifies a model variant which adds extra semantics of the
+[`ExpressionLimitedModel`](@ref), simulating the E-Flux algorithm. The
+arguments are forwarded to [`make_expression_limited_model`](@ref). Intended
+for usage with [`screen`](@ref).
+"""
+with_expression_limits(; kwargs...) =
+    model -> make_expression_limited_model(model; kwargs...)
diff --git a/test/analysis/expressions.jl b/test/analysis/expressions.jl
new file mode 100644
index 000000000..bbd5af8fc
--- /dev/null
+++ b/test/analysis/expressions.jl
@@ -0,0 +1,16 @@
+
+@testset "Expression-limited models" begin
+    orig_model = load_model(model_paths["e_coli_core.json"])
+
+    model =
+        orig_model |>
+        with_expression_limits(relative_expression = Dict(genes(orig_model) .=> 0.5))
+
+    bs = Dict(reactions(model) .=> bounds(model)[1])
+
+    @test getindex.(Ref(bs), ["PFK", "ENO", "ACALD"]) == [0, -500, -750]
+
+    fluxes = flux_balance_analysis_dict(model, Tulip.Optimizer)
+
+    @test isapprox(fluxes["BIOMASS_Ecoli_core_w_GAM"], 0.21386, atol = TEST_TOLERANCE)
+end

From eed8be5969eb56e7d9d327b27cca0e986d72dcf7 Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <miroslav.kratochvil@uni.lu>
Date: Wed, 3 May 2023 10:47:06 +0200
Subject: [PATCH 20/23] bump version

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 9f5b2465a..1d5443614 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "COBREXA"
 uuid = "babc4406-5200-4a30-9033-bf5ae714c842"
 authors = ["The developers of COBREXA.jl"]
-version = "1.4.4"
+version = "1.5"
 
 [deps]
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"

From 63200ef4a223277be9d9bc069c81a0a74b81a4cf Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <miroslav.kratochvil@uni.lu>
Date: Mon, 8 May 2023 09:30:27 +0200
Subject: [PATCH 21/23] fix filename

---
 src/reconstruction/{expresisons.jl => expressions.jl} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename src/reconstruction/{expresisons.jl => expressions.jl} (100%)

diff --git a/src/reconstruction/expresisons.jl b/src/reconstruction/expressions.jl
similarity index 100%
rename from src/reconstruction/expresisons.jl
rename to src/reconstruction/expressions.jl

From 391b84053c38b7e31c163587804154288a077e4c Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <miroslav.kratochvil@uni.lu>
Date: Thu, 11 May 2023 09:44:20 +0200
Subject: [PATCH 22/23] fix the change of default behavior of sparsevec() in
 1.9

---
 src/base/types/SBMLModel.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/base/types/SBMLModel.jl b/src/base/types/SBMLModel.jl
index 96c0633b8..afd13ebb0 100644
--- a/src/base/types/SBMLModel.jl
+++ b/src/base/types/SBMLModel.jl
@@ -180,7 +180,7 @@ $(TYPEDSIGNATURES)
 Objective of the [`SBMLModel`](@ref).
 """
 function objective(model::SBMLModel)::SparseVec
-    res = sparsevec([], [], n_reactions(model))
+    res = spzeros(n_reactions(model))
 
     objective = get(model.sbml.objectives, model.active_objective, nothing)
     if isnothing(objective) && length(model.sbml.objectives) == 1

From 1bebba915f604b43940d1e31528d3f546e0f4a3f Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <miroslav.kratochvil@uni.lu>
Date: Thu, 11 May 2023 09:45:19 +0200
Subject: [PATCH 23/23] release immediately

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 1d5443614..ee3a4ea68 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "COBREXA"
 uuid = "babc4406-5200-4a30-9033-bf5ae714c842"
 authors = ["The developers of COBREXA.jl"]
-version = "1.5"
+version = "1.5.1"
 
 [deps]
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"