Implement quadratic functionality

Project.toml

@@ -17,10 +17,10 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 
 [compat]
-AbstractFBCModels = "0.2"
+AbstractFBCModels = "0.2.2"
 Aqua = "0.7"
-Clarabel = "0.3"
-ConstraintTrees = "0.5"
+Clarabel = "0.6"
+ConstraintTrees = "0.6.1"
 DistributedData = "0.2"
 DocStringExtensions = "0.8, 0.9"
 Downloads = "1"
@@ -46,4 +46,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"
 
 [targets]
-test = ["Aqua", "Clarabel", "Downloads", "GLPK", "SHA", "Test", "Tulip", "JSONFBCModels"]
+test = ["Aqua", "Clarabel", "Downloads", "GLPK", "JSONFBCModels", "SBMLFBCModels", "SHA", "Test", "Tulip"]

docs/make.jl

@@ -2,6 +2,10 @@ using Documenter
 using Literate, JSON
 using COBREXA
 
+# testing constants
+const TEST_TOLERANCE = 1e-3
+const QP_TEST_TOLERANCE = 1e-2 # for Clarabel
+
 # build the examples
 examples_path = joinpath(@__DIR__, "src", "examples")
 examples_basenames = sort(filter(x -> endswith(x, ".jl"), readdir(examples_path)))

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

@@ -5,8 +5,6 @@
 # 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
-
 download_model(
     "http://bigg.ucsd.edu/static/models/e_coli_core.json",
     "e_coli_core.json",

docs/src/examples/02c-constraint-modifications.jl

@@ -43,7 +43,12 @@ fermentation = ctmodel.fluxes.EX_ac_e.value + ctmodel.fluxes.EX_etoh_e.value
 forced_mixed_fermentation =
     ctmodel * :fermentation^C.Constraint(fermentation, (10.0, 1000.0)) # new modified model is created
 
-vt = flux_balance(forced_mixed_fermentation, Tulip.Optimizer; modifications = [silence])
+vt = optimized_constraints(
+    forced_mixed_fermentation,
+    objective = forced_mixed_fermentation.objective.value,
+    optimizer = Tulip.Optimizer,
+    modifications = [silence],
+)
 
 @test isapprox(vt.objective, 0.6337, atol = TEST_TOLERANCE) #src
 
@@ -52,13 +57,25 @@ vt = flux_balance(forced_mixed_fermentation, Tulip.Optimizer; modifications = [s
 
 ctmodel.fluxes.ATPM.bound = (1000.0, 10000.0)
 
-vt = flux_balance(ctmodel, Tulip.Optimizer; modifications = [silence])
+#TODO explicitly show here how false sharing looks like
+
+vt = optimized_constraints(
+    ctmodel,
+    objective = ctmodel.objective.value,
+    optimizer = Tulip.Optimizer,
+    modifications = [silence],
+)
 
 @test isnothing(vt) #src
 
 # Models can also be piped into the analysis functions
 
 ctmodel.fluxes.ATPM.bound = (8.39, 10000.0) # revert
-vt = ctmodel |> flux_balance(Tulip.Optimizer; modifications = [silence])
+vt = optimized_constraints(
+    ctmodel,
+    objective = ctmodel.objective.value,
+    optimizer = Tulip.Optimizer,
+    modifications = [silence],
+)
 
 @test isapprox(vt.objective, 0.8739, atol = TEST_TOLERANCE) #src

docs/src/examples/03-parsimonious-flux-balance.jl

@@ -0,0 +1,99 @@
+
+# # Parsimonious flux balance analysis
+
+# We will use [`parsimonious_flux_balance`](@ref) and
+# [`minimize_metabolic_adjustment`](@ref) to find the optimal flux
+# distribution in the *E. coli* "core" model.
+#
+# TODO pFBA citation
+
+# If it is not already present, download the model and load the package:
+import Downloads: download
+
+!isfile("e_coli_core.json") &&
+    download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json")
+
+# next, load the necessary packages
+
+using COBREXA
+
+import JSONFBCModels
+import Clarabel # can solve QPs
+
+model = load_model("e_coli_core.json") # load the model
+
+# Use the convenience function to run standard pFBA on
+
+vt = parsimonious_flux_balance(model, Clarabel.Optimizer; modifications = [silence])
+
+# Or use the piping functionality
+
+model |> parsimonious_flux_balance(Clarabel.Optimizer; modifications = [silence])
+
+@test isapprox(vt.objective, 0.87392; atol = TEST_TOLERANCE) #src
+@test sum(x^2 for x in values(vt.fluxes)) < 15000 #src
+
+#=
+
+# Alternatively, you can construct your own constraint tree model with
+# the quadratic objective (this approach is much more flexible).
+
+ctmodel = fbc_model_constraints(model)
+ctmodel *= :l2objective^squared_sum_objective(ctmodel.fluxes)
+ctmodel.objective.bound = 0.3 # set growth rate # TODO currently breaks
+
+opt_model = optimization_model(
+    ctmodel;
+    objective = ctmodel.:l2objective.value,
+    optimizer = Clarabel.Optimizer,
+    sense = Minimal,
+)
+
+J.optimize!(opt_model) # JuMP is called J in COBREXA
+
+is_solved(opt_model) # check if solved
+
+vt = C.constraint_values(ctmodel, J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA
+
+@test isapprox(vt.l2objective, ?; atol = QP_TEST_TOLERANCE) #src  # TODO will break until mutable bounds
+
+# It is likewise as simple to run MOMA using the convenience functions.
+
+ref_sol = Dict("ATPS4r" => 33.0, "CYTBD" => 22.0)
+
+vt = minimize_metabolic_adjustment(model, ref_sol, Gurobi.Optimizer)
+
+# Or use the piping functionality
+
+model |>
+minimize_metabolic_adjustment(ref_sol, Clarabel.Optimizer; modifications = [silence])
+
+@test isapprox(vt.:momaobjective, 0.81580806; atol = TEST_TOLERANCE) #src
+
+# Alternatively, you can construct your own constraint tree model with
+# the quadratic objective (this approach is much more flexible).
+
+ctmodel = fbc_model_constraints(model)
+ctmodel *=
+    :minoxphospho^squared_sum_error_objective(
+        ctmodel.fluxes,
+        Dict(:ATPS4r => 33.0, :CYTBD => 22.0),
+    )
+ctmodel.objective.bound = 0.3 # set growth rate # TODO currently breaks
+
+opt_model = optimization_model(
+    ctmodel;
+    objective = ctmodel.minoxphospho.value,
+    optimizer = Clarabel.Optimizer,
+    sense = Minimal,
+)
+
+J.optimize!(opt_model) # JuMP is called J in COBREXA
+
+is_solved(opt_model) # check if solved
+
+vt = C.constraint_values(ctmodel, J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA
+
+@test isapprox(vt.l2objective, ?; atol = QP_TEST_TOLERANCE) #src  # TODO will break until mutable bounds
+
+=#

docs/src/examples/04-minimization-of-metabolic-adjustment.jl

@@ -0,0 +1,23 @@
+
+# # Minimization of metabolic adjustment
+
+# TODO MOMA citation
+
+import Downloads: download
+
+!isfile("e_coli_core.json") &&
+    download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json")
+
+using COBREXA
+import AbstractFBCModels.CanonicalModel as CM
+import JSONFBCModels
+import Clarabel
+
+# TODO this might do the convert immediately as with the old cobrexa...
+# probably better have an actual output-type argument tho rather than force the
+# guessing.
+model = convert(CM.Model, load_model("e_coli_core.json"))
+
+reference_fluxes = parsimonious_flux_balance(model, Clarabel.Optimizer).fluxes
+
+# TODO MOMA from here

src/COBREXA.jl

@@ -8,6 +8,7 @@ using DocStringExtensions
 import AbstractFBCModels as A
 import ConstraintTrees as C
 import JuMP as J
+import SparseArrays: sparse
 
 include("types.jl")
 include("io.jl")
@@ -20,5 +21,8 @@ include("builders/objectives.jl")
 include("analysis/modifications.jl")
 include("analysis/flux_balance.jl")
 include("analysis/parsimonious_flux_balance.jl")
+include("analysis/minimal_metabolic_adjustment.jl")
+
+include("misc/bounds.jl")
 
 end # module COBREXA

src/analysis/flux_balance.jl

@@ -1,64 +1,54 @@
 """
 $(TYPEDSIGNATURES)
 
-Run flux balance analysis (FBA) on the `model`, optionally specifying
-`modifications` to the problem.  Basically, FBA solves this optimization
-problem:
-```
-max cᵀx
-s.t. S x = b
-     xₗ ≤ x ≤ xᵤ
-```
-See "Orth, J., Thiele, I. & Palsson, B. What is flux balance analysis?. Nat
-Biotechnol 28, 245-248 (2010). https://doi.org/10.1038/nbt.1614" for more
-information.
+Make an JuMP model out of `constraints` using [`optimization_model`](@ref)
+(most arguments are forwarded there), then apply the `modifications`, optimize
+the model, and return either `nothing` if the optimization failed, or `output`
+substituted with the solved values (`output` defaults to `constraints`.
 
-The `optimizer` must be set to a `JuMP`-compatible optimizer, such as
-`GLPK.Optimizer` or `Tulip.Optimizer`.
-
-Optionally, you may specify one or more modifications to be applied to the model
-before the analysis, such as [`set_objective_sense`](@ref),
-[`set_optimizer`](@ref), [`set_optimizer_attribute`](@ref), and
-[`silence`](@ref).
-
-Returns a [`C.ValueTree`](@ref).
-
-# Example
-```
-model = load_model("e_coli_core.json")
-solution = flux_balance(model, GLPK.optimizer)
-```
+For a "nice" version for simpler finding of metabolic model optima, use
+[`flux_balance`](@ref).
 """
-function flux_balance(model::A.AbstractFBCModel, optimizer; modifications = [])
-    ctmodel = fbc_model_constraints(model)
-    flux_balance(ctmodel, optimizer; modifications)
+function optimized_constraints(
+    constraints::C.ConstraintTreeElem;
+    modifications = [],
+    output = constraints,
+    kwargs...,
+)
+    om = optimization_model(constraints; kwargs...)
+    for m in modifications
+        m(om)
+    end
+    J.optimize!(om)
+    is_solved(om) ? C.constraint_values(output, J.value.(om[:x])) : nothing
 end
 
-"""
-$(TYPEDSIGNATURES)
+export optimized_constraints
 
-A variant of [`flux_balance`](@ref) that takes in a
-[`C.ConstraintTree`](@ref) as the model to optimize. The objective is inferred
-from the field `objective` in `ctmodel`. All other arguments are forwarded.
 """
-function flux_balance(ctmodel::C.ConstraintTree, optimizer; modifications = [])
-    opt_model = optimization_model(ctmodel; objective = ctmodel.objective.value, optimizer)
-
-    for mod in modifications
-        mod(opt_model)
-    end
+$(TYPEDSIGNATURES)
 
-    J.optimize!(opt_model)
+Compute an optimal objective-optimizing solution of the given `model`.
 
-    is_solved(opt_model) || return nothing
+Most arguments are forwarded to [`optimized_constraints`](@ref).
 
-    C.ValueTree(ctmodel, J.value.(opt_model[:x]))
+Returns a tree with the optimization solution of the same shape as
+given by [`fbc_model_constraints`](@ref).
+"""
+function flux_balance(model::A.AbstractFBCModel, optimizer; kwargs...)
+    constraints = fbc_model_constraints(model)
+    optimized_constraints(
+        constraints;
+        objective = constraints.objective.value,
+        optimizer,
+        kwargs...,
+    )
 end
 
 """
 $(TYPEDSIGNATURES)
 
-Pipe-able variant of [`flux_balance`](@ref).
+Pipe-able overload of [`flux_balance`](@ref).
 """
 flux_balance(optimizer; modifications = []) = m -> flux_balance(m, optimizer; modifications)
 

src/analysis/minimal_metabolic_adjustment.jl

@@ -0,0 +1,77 @@
+
+"""
+$(TYPEDSIGNATURES)
+
+Find a feasible solution of the "minimal metabolic adjustment analysis" (MOMA)
+for the `model`, which is the "closest" feasible solution to the given
+`reference_fluxes`, in the sense of squared-sum error distance. The minimized
+squared distance (the objective) is present in the result tree as
+`minimal_adjustment_objective`.
+
+This is often used for models with smaller feasible region than the reference
+models (typically handicapped by a knockout, nutritional deficiency or a
+similar perturbation). MOMA solution then gives an expectable "easiest"
+adjustment of the organism towards a somewhat working state.
+
+Reference fluxes that do not exist in the model are ignored (internally, the
+objective is constructed via [`squared_sum_error_objective`](@ref)).
+
+Additional parameters are forwarded to [`optimized_constraints`](@ref).
+"""
+function minimal_metabolic_adjustment(
+    model::A.AbstractFBCModel,
+    reference_fluxes::Dict{Symbol,Float64},
+    optimizer;
+    kwargs...,
+)
+    constraints = fbc_model_constraints(model)
+    objective = squared_sum_error_objective(constraints.fluxes, reference_fluxes)
+    optimized_constraints(
+        constraints * :minimal_adjustment_objective^C.Constraint(objective);
+        optimizer,
+        objective,
+        sense = Minimal,
+        kwargs...,
+    )
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+A slightly easier-to-use version of [`minimal_metabolic_adjustment`](@ref) that
+computes the reference flux as the optimal solution of the
+[`reference_model`](@ref). The reference flux is calculated using
+`reference_optimizer` and `reference_modifications`, which default to the
+`optimizer` and `modifications`.
+
+Leftover arguments are passed to the overload of
+[`minimal_metabolic_adjustment`](@ref) that accepts the reference flux
+dictionary.
+"""
+function minimal_metabolic_adjustment(
+    model::A.AbstractFBCModel,
+    reference_model::A.AbstractFBCModel,
+    optimizer;
+    reference_optimizer = optimizer,
+    modifications = [],
+    reference_modifications = modifications,
+    kwargs...,
+)
+    reference_constraints = fbc_model_constraints(reference_model)
+    reference_fluxes = optimized_constraints(
+        reference_constraints;
+        optimizer = reference_optimizer,
+        modifications = reference_modifications,
+        output = reference_constraints.fluxes,
+    )
+    isnothing(reference_fluxes) && return nothing
+    minimal_metabolic_adjustment(
+        model,
+        reference_fluxes,
+        optimizer;
+        modifications,
+        kwargs...,
+    )
+end
+
+export minimal_metabolic_adjustment