Merge pull request #800 from LCSB-BioCore/sew-fba-tests-etc

Start getting FBA to work

Author: exaexa

.gitignore

@@ -20,3 +20,8 @@ Manifest.toml
 
 # ignore container files
 *.sif
+
+# ignore models
+*.xml
+*.json
+*.mat

.gitlab-ci.yml

@@ -71,9 +71,9 @@ variables:
   before_script:
     - docker login -u $CI_USER_NAME -p $GITLAB_ACCESS_TOKEN $CI_REGISTRY
 
-.global_julia18: &global_julia18
+.global_julia19: &global_julia19
   variables:
-    JULIA_VER: "v1.8.3"
+    JULIA_VER: "v1.9.4"
 
 .global_julia16: &global_julia16
   variables:
@@ -120,26 +120,24 @@ variables:
 # any available docker and current julia
 #
 
-docker:julia1.8:
+docker:julia1.9:
   stage: test
   image: $CI_REGISTRY/r3/docker/julia-custom
   script:
     - julia --check-bounds=yes --inline=yes --project=@. -e "import Pkg; Pkg.test(; coverage = true)"
-  after_script:
-    - julia --project=test/coverage test/coverage/coverage-summary.jl
   <<: *global_trigger_pull_request
 
 #
 # The required compatibility test to pass on branches&tags before the docs get
 # built & deployed
 #
 
-linux:julia1.8:
+linux:julia1.9:
   stage: test
   tags:
     - slave01
   <<: *global_trigger_full_tests
-  <<: *global_julia18
+  <<: *global_julia19
   <<: *global_env_linux
 
 linux:julia1.6:
@@ -154,22 +152,22 @@ linux:julia1.6:
 # Additional platform&environment compatibility tests
 #
 
-windows8:julia1.8:
+windows8:julia1.9:
   stage: test-compat
   <<: *global_trigger_compat_tests
-  <<: *global_julia18
+  <<: *global_julia19
   <<: *global_env_win8
 
-windows10:julia1.8:
+windows10:julia1.9:
   stage: test-compat
   <<: *global_trigger_compat_tests
-  <<: *global_julia18
+  <<: *global_julia19
   <<: *global_env_win10
 
-mac:julia1.8:
+mac:julia1.9:
   stage: test-compat
   <<: *global_trigger_compat_tests
-  <<: *global_julia18
+  <<: *global_julia19
   <<: *global_env_mac
 
 windows8:julia1.6:

Project.toml

@@ -18,20 +18,32 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 
 [compat]
 AbstractFBCModels = "0.2"
+Aqua = "0.7"
 Clarabel = "0.3"
-ConstraintTrees = "0.4"
+ConstraintTrees = "0.5"
 DistributedData = "0.2"
 DocStringExtensions = "0.8, 0.9"
+Downloads = "1"
+GLPK = "1"
+JSONFBCModels = "0.1"
 JuMP = "1"
+SBMLFBCModels = "0.1"
+SHA = "0.7, 1"
 StableRNGs = "1.0"
+Test = "1"
+Tulip = "0.9"
 julia = "1.5"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
+Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
+JSONFBCModels = "475c1105-d6ed-49c1-9b32-c11adca6d3e8"
+SBMLFBCModels = "3e8f9d1a-ffc1-486d-82d6-6c7276635980"
+SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"
 
 [targets]
-test = ["Aqua", "Clarabel", "GLPK", "Test", "Tulip"]
+test = ["Aqua", "Clarabel", "Downloads", "GLPK", "SHA", "Test", "Tulip", "JSONFBCModels"]

README.md

@@ -115,7 +115,7 @@ download("http://bigg.ucsd.edu/static/models/e_coli_core.xml", "e_coli_core.xml"
 model = load_model("e_coli_core.xml")
 
 # run a FBA
-fluxes = flux_balance_analysis_dict(model, Tulip.Optimizer)
+fluxes = flux_balance_dict(model, Tulip.Optimizer)
 ```
 
 The variable `fluxes` will now contain a dictionary of the computed optimal

docs/Project.toml

@@ -1,15 +1,10 @@
 [deps]
+AbstractFBCModels = "5a4f3dfa-1789-40f8-8221-69268c29937c"
 COBREXA = "babc4406-5200-4a30-9033-bf5ae714c842"
-CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
-Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
-ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-Escher = "8cc96de1-1b23-48cb-9272-618d67962629"
-GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+JSONFBCModels = "475c1105-d6ed-49c1-9b32-c11adca6d3e8"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"
 

docs/make.jl

@@ -25,6 +25,7 @@ find_mds(path) =
         filter(x -> endswith(x, ".md"), readdir(joinpath(@__DIR__, "src", path))),
     )
 
+#TODO migrate this to Documenter-1, and make all checks strict
 # build the docs
 makedocs(
     modules = [COBREXA],
@@ -51,10 +52,11 @@ makedocs(
             "Contents" => "concepts.md"
             find_mds("concepts")
         ],
-        "Reference" => [
-            "Contents" => "reference.md"
-            find_mds("reference")
-        ],
+        "Reference" => "reference.md",
+        #[ # TODO re-add this when the reference gets bigger
+        #"Contents" => "reference.md"
+        #find_mds("reference")
+        #],
     ],
 )
 

docs/src/concepts/1_screen.md

@@ -24,7 +24,7 @@ screen_variants(
         [with_changed_bound("O2t", lb = 0, ub = 0)],  # disable O2 transport
         [with_changed_bound("CO2t", lb = 0, ub = 0), with_changed_bound("O2t", lb = 0, ub = 0)],  # disable both transports
     ],
-    m -> flux_balance_analysis_dict(m, Tulip.Optimizer)["BIOMASS_Ecoli_core_w_GAM"],
+    m -> flux_balance_dict(m, Tulip.Optimizer)["BIOMASS_Ecoli_core_w_GAM"],
 )
 ```
 The call specifies a model (the `m` that we have loaded) that is being tested,
@@ -89,7 +89,7 @@ res = screen_variants(m,
             ["EX_h2o_e", "EX_co2_e", "EX_o2_e", "EX_nh4_e"], # and this set of exchanges
         )
     ],
-    m -> flux_balance_analysis_dict(m, Tulip.Optimizer)["BIOMASS_Ecoli_core_w_GAM"],
+    m -> flux_balance_dict(m, Tulip.Optimizer)["BIOMASS_Ecoli_core_w_GAM"],
 )
 ```
 
@@ -199,7 +199,7 @@ screen_variants(
         [with_disabled_oxygen_transport],
         [with_disabled_reaction("NH4t")],
     ],
-    m -> flux_balance_analysis_dict(m, Tulip.Optimizer)["BIOMASS_Ecoli_core_w_GAM"],
+    m -> flux_balance_dict(m, Tulip.Optimizer)["BIOMASS_Ecoli_core_w_GAM"],
 )
 ```
 
@@ -222,7 +222,7 @@ That should get you the results for all new variants of the model:
 
 Some analysis functions may take additional arguments, which you might want to
 vary for the analysis. `modifications` argument of
-[`flux_balance_analysis_dict`](@ref) is one example of such argument, allowing
+[`flux_balance_dict`](@ref) is one example of such argument, allowing
 you to specify details of the optimization procedure.
 
 [`screen`](@ref) function allows you to do precisely that -- apart from
@@ -242,7 +242,7 @@ iterations needed for Tulip solver to find a feasible solution:
 screen(m,
     args = [(i,) for i in 5:15],  # the iteration counts, packed in 1-tuples
     analysis = (m,a) -> # `args` elements get passed as the extra parameter here
-        flux_balance_analysis_vec(m,
+        flux_balance_vec(m,
             Tulip.Optimizer;
             modifications=[modify_optimizer_attribute("IPM_IterationsLimit", a)],
         ),

docs/src/examples/01-loading-and-saving.jl

@@ -4,3 +4,5 @@
 using COBREXA
 
 # TODO: download the models into a single directory that can get cached. Probably best have a fake mktempdir().
+#
+# TODO: demonstrate download_model here and explain how to get hashes (simply not fill them in for the first time)

docs/src/examples/02-flux-balance-analysis.jl

@@ -1,6 +1,51 @@
 
-# # Flux balance analysis
+# # Flux balance analysis (FBA)
+
+# Here we use [`flux_balance`](@ref) and several related functions to
+# find an optimal flux in the *E. coli* "core" model. We will need the model,
+# which we can download using [`download_model`](@ref):
 
 using COBREXA
 
-# TODO: run FBA on a FBC model
+download_model(
+    "http://bigg.ucsd.edu/static/models/e_coli_core.json",
+    "e_coli_core.json",
+    "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
+)
+
+# Additionally to COBREXA and the model format package, we will need a solver
+# -- let's use Tulip here:
+
+import JSONFBCModels
+import Tulip
+
+model = load_model("e_coli_core.json")
+
+# ## Running a FBA
+#
+# There are many possibilities on how to arrange the metabolic model into the
+# optimization framework and how to actually solve it. The "usual" assumed one
+# is captured in the default behavior of function
+# [`flux_balance`](@ref):
+
+solution = flux_balance(model, Tulip.Optimizer)
+
+@test isapprox(solution.objective, 0.8739, atol = TEST_TOLERANCE) #src
+
+# The result contains a tree of all optimized values in the model, including
+# fluxes, the objective value, and possibly others (given by what the model
+# contains).
+#
+# You can explore the dot notation to explore the solution, extracting e.g. the
+# value of the objective:
+
+solution.objective
+
+# ...or the value of the flux through the given reaction (note the solution is
+# not unique in FBA):
+
+solution.fluxes.PFK
+
+# ...or make a "table" of all fluxes through all reactions:
+
+collect(solution.fluxes)

docs/src/examples/02a-optimizer-parameters.jl

@@ -0,0 +1,57 @@
+
+# # Changing optimizer parameters
+#
+# Many optimizers require fine-tuning to produce best results. You can pass in
+# additional optimizer settings via the `modifications` parameter of
+# [`flux_balance`](@ref). These include e.g.
+#
+# - [`set_optimizer_attribute`](@ref) (typically allowing you to tune e.g.
+#   iteration limits, tolerances, or floating-point precision)
+# - [`set_objective_sense`](@ref) (allowing you to change and reverse the
+#   optimization direction, if required)
+# - [`silence`](@ref) to disable the debug output of the optimizer
+# - and even [`set_optimizer`](@ref), which changes the optimizer
+#   implementation used (this is not quite useful in this case, but becomes
+#   beneficial with more complex, multi-stage optimization problems)
+#
+# To demonstrate this, we'll use the usual toy model:
+
+using COBREXA
+import JSONFBCModels, Tulip
+
+download_model(
+    "http://bigg.ucsd.edu/static/models/e_coli_core.json",
+    "e_coli_core.json",
+    "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
+)
+
+model = load_model("e_coli_core.json")
+
+# Running a FBA with a silent optimizer that has slightly increased iteration
+# limit for IPM algorithm may now look as follows:
+solution = flux_balance(
+    model,
+    Tulip.Optimizer;
+    modifications = [silence, set_optimizer_attribute("IPM_IterationsLimit", 1000)],
+)
+
+@test !isnothing(solution) #src
+
+# To see some of the effects of the configuration changes, you may e.g.
+# deliberately cripple the optimizer's possibilities to a few iterations, which
+# will cause it to fail and return no solution:
+
+solution = flux_balance(
+    model,
+    Tulip.Optimizer;
+    modifications = [silence, set_optimizer_attribute("IPM_IterationsLimit", 2)],
+)
+
+println(solution)
+
+@test isnothing(solution) #src
+
+# Applicable optimizer attributes are documented in the documentations of the
+# respective optimizers. To browse the possibilities, you may want to see the
+# [JuMP documentation page that summarizes the references to the available
+# optimizers](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers).