Skip to content

Commit

Permalink
fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
exaexa committed Dec 9, 2023
1 parent 48cbc6e commit a2c121e
Show file tree
Hide file tree
Showing 4 changed files with 58 additions and 29 deletions.
6 changes: 3 additions & 3 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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.6"
Clarabel = "0.6"
ConstraintTrees = "0.6.1"
DistributedData = "0.2"
DocStringExtensions = "0.8, 0.9"
Downloads = "1"
Expand Down
Original file line number Diff line number Diff line change
@@ -1,9 +1,11 @@

# # Quandratic objective flux balance analysis type problems
# # 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
Expand All @@ -13,81 +15,85 @@ import Downloads: download

# next, load the necessary packages

import COBREXA as X
import AbstractFBCModels as A # for the accessors
import JSONFBCModels as J # for the model type
using COBREXA

import JSONFBCModels
import Clarabel # can solve QPs

model = A.load(J.JSONFBCModel, "e_coli_core.json") # load the model
model = load_model("e_coli_core.json") # load the model

# Use the convenience function to run standard pFBA
# Use the convenience function to run standard pFBA on

vt = X.parsimonious_flux_balance(model, Clarabel.Optimizer; modifications = [X.silence])
vt = parsimonious_flux_balance(model, Clarabel.Optimizer; modifications = [silence])

# Or use the piping functionality

model |> parsimonious_flux_balance(Clarabel.Optimizer; modifications = [X.silence])
model |> parsimonious_flux_balance(Clarabel.Optimizer; modifications = [silence])

@test isapprox(vt.objective, 0.87392; atol = TEST_TOLERANCE) #src
@test isapprox(sum(x^2 for x in values(vt.fluxes)), 11414.2119; atol = QP_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 = X.fbc_model_constraints(model)
ctmodel *= :l2objective^X.squared_sum_objective(ctmodel.fluxes)
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 = X.optimization_model(
opt_model = optimization_model(
ctmodel;
objective = ctmodel.:l2objective.value,
optimizer = Clarabel.Optimizer,
sense = X.J.MIN_SENSE,
sense = J.MIN_SENSE,
)
X.J.optimize!(opt_model) # JuMP is called J in COBREXA
J.optimize!(opt_model) # JuMP is called J in COBREXA
X.is_solved(opt_model) # check if solved
is_solved(opt_model) # check if solved
vt = X.C.constraint_values(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA
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 = X.minimize_metabolic_adjustment(model, ref_sol, Gurobi.Optimizer)
vt = minimize_metabolic_adjustment(model, ref_sol, Gurobi.Optimizer)
# Or use the piping functionality
model |>
X.minimize_metabolic_adjustment(ref_sol, Clarabel.Optimizer; modifications = [X.silence])
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 = X.fbc_model_constraints(model)
ctmodel = fbc_model_constraints(model)
ctmodel *=
:minoxphospho^X.squared_sum_error_objective(
: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 = X.optimization_model(
opt_model = optimization_model(
ctmodel;
objective = ctmodel.minoxphospho.value,
optimizer = Clarabel.Optimizer,
sense = X.J.MIN_SENSE,
sense = J.MIN_SENSE,
)
X.J.optimize!(opt_model) # JuMP is called J in COBREXA
J.optimize!(opt_model) # JuMP is called J in COBREXA
X.is_solved(opt_model) # check if solved
is_solved(opt_model) # check if solved
vt = X.C.constraint_values(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA
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
=#
23 changes: 23 additions & 0 deletions docs/src/examples/04-minimization-of-metabolic-adjustment.jl
Original file line number Diff line number Diff line change
@@ -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
2 changes: 1 addition & 1 deletion src/analysis/parsimonious_flux_balance.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ function parsimonious_optimized_constraints(
)

# first solve the optimization problem with the original objective
om = optimization_model(constraints, args...; kwargs...)
om = optimization_model(constraints, args...; objective, kwargs...)
for m in modifications
m(om)
end
Expand Down

0 comments on commit a2c121e

Please sign in to comment.