diff --git a/src/analysis/fba.jl b/src/analysis/fba.jl
index 2e15bf838e2017034e8883ab35aa34d1df4d427a..f5adcfdd363e57eeed2b5bfa0964bb7f9939f4e7 100644
--- a/src/analysis/fba.jl
+++ b/src/analysis/fba.jl
@@ -21,7 +21,8 @@ of the model, if the solution is found.
Arguments are passed to [`flux_balance_analysis`](@ref).
"""
function flux_balance_analysis_vec(args...)::Union{Vector{Float64},Nothing}
- (optmodel, vars) = flux_balance_analysis(args...)
+ optmodel = flux_balance_analysis(args...)
+ vars = optmodel[:x]
termination_status(optmodel) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
value.(vars)
@@ -71,7 +72,8 @@ function fba(
sense = MOI.MAX_SENSE,
)
# get core optimization problem
- cbm, v, mb, lbcons, ubcons = make_optimization_model(model, optimizer, sense = sense)
+ cbm = make_optimization_model(model, optimizer, sense = sense)
+ v = cbm[:x] # fluxes
# modify core optimization problem according to user specifications
if !isempty(solver_attributes) # set other attributes
@@ -83,7 +85,7 @@ function fba(
# set additional constraints
for (rxnid, con) in constraints
ind = model.reactions[findfirst(model.reactions, rxnid)]
- set_bound(ind, lbcons, ubcons; lb = con[1], ub = con[2])
+ set_bound(ind, cbm; lb = con[1], ub = con[2])
end
# if an objective function is supplied, modify the default objective
diff --git a/src/analysis/fva.jl b/src/analysis/fva.jl
index 0b38dc87179e3d82db28c96b2f5abcd4db8a0f6c..e1c5a575c1470b3341423f9d3095d1cd532efda9 100644
--- a/src/analysis/fva.jl
+++ b/src/analysis/fva.jl
@@ -37,14 +37,15 @@ function flux_variability_analysis(
throw(DomainError(reactions, "Index exceeds number of reactions."))
end
- (optimization_model, x0) = flux_balance_analysis(model, optimizer)
+ optimization_model = flux_balance_analysis(model, optimizer)
Z0 = COBREXA.JuMP.objective_value(optimization_model)
optimization_model = nothing # we won't need this one anymore, so free the memory
# store a JuMP optimization model at all workers
save_model = :(
begin
- optmodel, x = COBREXA.make_optimization_model($model, $optimizer)
+ optmodel = COBREXA.make_optimization_model($model, $optimizer)
+ x = optmodel[:x]
COBREXA._FVA_add_constraint(optmodel, $(objective(model)), x, $Z0, $gamma)
optmodel
end
@@ -140,7 +141,8 @@ function fva(
constraints = Dict{String,Tuple{Float64,Float64}}(),
sense = MOI.MAX_SENSE,
)
- cbm, v, mb, lbcons, ubcons = make_optimization_model(model, optimizer, sense = sense)
+ cbm = make_optimization_model(model, optimizer, sense = sense)
+ v = cbm[:x]
if !isempty(solver_attributes) # set other attributes
for (k, v) in solver_attributes
@@ -151,7 +153,7 @@ function fva(
# set additional constraints
for (rxnid, con) in constraints
ind = model.reactions[findfirst(model.reactions, rxnid)]
- set_bound(ind, lbcons, ubcons; lb = con[1], ub = con[2])
+ set_bound(ind, cbm; lb = con[1], ub = con[2])
end
# if an objective function is supplied, modify the default objective
diff --git a/src/analysis/pfba.jl b/src/analysis/pfba.jl
index 1e2afc4139a4494610cda97493d8da0643a5a941..8402d3f4da2052a72ddc9f9003a211a1d44c4047 100644
--- a/src/analysis/pfba.jl
+++ b/src/analysis/pfba.jl
@@ -33,16 +33,16 @@ function pfba(
if typeof(optimizer) <: AbstractArray # choose optimizer
- cbm, v, mb, lbcons, ubcons =
- make_optimization_model(model, optimizer[1], sense = sense)
+ cbm = make_optimization_model(model, optimizer[1], sense = sense)
+ v = cbm[:x]
if !isempty(solver_attributes["opt1"]) # set other attributes
for (k, v) in solver_attributes["opt1"]
set_optimizer_attribute(cbm, k, v)
end
end
else # singe optimizer
- cbm, v, mb, lbcons, ubcons =
- make_optimization_model(model, optimizer, sense = sense)
+ cbm = make_optimization_model(model, optimizer, sense = sense)
+ v = cbm[:x]
if !isempty(solver_attributes) # set other attributes
for (k, v) in solver_attributes
set_optimizer_attribute(cbm, k, v)
@@ -53,7 +53,7 @@ function pfba(
# set additional constraints
for (rxnid, con) in constraints
ind = model.reactions[findfirst(model.reactions, rxnid)]
- set_bound(ind, lbcons, ubcons; lb = con[1], ub = con[2])
+ set_bound(ind, cbm; lb = con[1], ub = con[2])
end
# check if default objective should be used
diff --git a/src/base/solver.jl b/src/base/solver.jl
index 178654463f73c57a05844438ca73c0cdf8674620..b8d16887170516b01d66a9bc8d4e31c08e5f8a66 100644
--- a/src/base/solver.jl
+++ b/src/base/solver.jl
@@ -13,11 +13,11 @@ function make_optimization_model(
optimization_model = COBREXA.JuMP.Model(optimizer)
@variable(optimization_model, x[i = 1:n])
@objective(optimization_model, sense, objective(model)' * x)
- mb = @constraint(optimization_model, mb, stoichiometry(model) * x .== balance(model)) # mass balance
- lbs = @constraint(optimization_model, lbs, xl .<= x) # lower bounds
- ubs = @constraint(optimization_model, ubs, x .<= xu) # upper bounds
+ @constraint(optimization_model, mb, stoichiometry(model) * x .== balance(model)) # mass balance
+ @constraint(optimization_model, lbs, xl .<= x) # lower bounds
+ @constraint(optimization_model, ubs, x .<= xu) # upper bounds
- return optimization_model, x, mb, lbs, ubs
+ return optimization_model
end
"""
@@ -28,8 +28,7 @@ function optimize_model(
optimizer;
sense = MOI.MIN_SENSE,
) where {LM<:MetabolicModel}
- optimization_model, x, _, _, _ =
- make_optimization_model(model, optimizer; sense = sense)
- COBREXA.JuMP.optimize!(optimization_model)
- return (optimization_model, x)
+ optimization_model = make_optimization_model(model, optimizer; sense = sense)
+ COBREXA.JuMP.optimize!(optimization_model)
+ return optimization_model
end
diff --git a/src/base/utilities.jl b/src/base/utilities.jl
index 6ea40a1767aae5e13e9b23f6a48d0c07c804d242..ac473dd5677cb0e740583539a7939ae4fabd251d 100644
--- a/src/base/utilities.jl
+++ b/src/base/utilities.jl
@@ -21,16 +21,16 @@ Base.copy(model::CoupledLinearModel) =
"""
- set_bound(index, ubconstaintref, lbconstaintref; ub=1000, lb=-1000)
+ set_bound(index, optimization_model; ub=1000, lb=-1000)
Helper function to set the bounds of variables.
The JuMP `set_normalized_rhs` function is a little confusing, so this function simplifies setting
-constraints. Just supply the constraint `index` and the bound objects (`ubconstaintref`, `lbconstaintref`) and they will be set to `ub` and `lb`.
+constraints. Just supply the constraint `index` and the model and they will be set to `ub` and `lb`.
"""
-function set_bound(vind, lbs, ubs; ub = 1000, lb = -1000)
+function set_bound(vind, opt_model; ub = 1000, lb = -1000)
if lb <= 0
- set_normalized_rhs(lbs[vind], abs(lb))
+ set_normalized_rhs(opt_model[:lbs][vind], abs(lb))
else
- set_normalized_rhs(lbs[vind], -abs(lb))
+ set_normalized_rhs(opt_model[:lbs][vind], -abs(lb))
end
- set_normalized_rhs(ubs[vind], ub)
+ set_normalized_rhs(opt_model[:ubs][vind], ub)
end
diff --git a/src/sampling/sampling_tools.jl b/src/sampling/sampling_tools.jl
index 0762c4cf5af0d6560963b69cd7542417fbe6c202..083c53dcd967389113a44e3cca42ae9c04a26371 100644
--- a/src/sampling/sampling_tools.jl
+++ b/src/sampling/sampling_tools.jl
@@ -1,14 +1,17 @@
"""
- get_warmup_points(cbmodel, v, mb, lbs, ubs; random_objective=false, numstop=1e10)
+ get_warmup_points(cbm; random_objective=false, numstop=1e10)
-Generate warmup points for all the reactions on the model that
+Generate warmup points for all the reactions in the model that
are not fixed. Assumes you feed in a JuMP model that is already
-constrained i.e. the constrains are already applied into cbmodel.
-Note, extra constraints applied to ubs and lbs will have no effect.
+constrained i.e. the constrains are already applied into `cbm`.
-numstop = 2*number of warmup points - to reduce the time this takes
+numstop (maximum value is 2*number of warmup points). Reduce this to prematurely
+stop finding warmup points.
"""
-function get_warmup_points(cbmodel, v, lbs, ubs; random_objective = false, numstop = 1e10)
+function get_warmup_points(cbm; random_objective = false, numstop = 1e10)
+ v = cbm[:x]
+ ubs = cbm[:ubs]
+ lbs = cbm[:lbs]
# determine which rxns should be max/min-ized
fixed_rxns = Int64[]
for i in eachindex(v)
@@ -30,23 +33,23 @@ function get_warmup_points(cbmodel, v, lbs, ubs; random_objective = false, numst
i > NN && break
if random_objective
- @objective(cbmodel, Max, sum(rand() * v[iii] for iii in var_rxn_inds))
+ @objective(cbm, Max, sum(rand() * v[iii] for iii in var_rxn_inds))
else
- @objective(cbmodel, Max, v[var_rxn_inds[i]])
+ @objective(cbm, Max, v[var_rxn_inds[i]])
end
- optimize!(cbmodel)
+ optimize!(cbm)
for j = 1:size(wpoints, 1)
wpoints[j, ii] = value(v[j])
end
if random_objective
- @objective(cbmodel, Min, sum(rand() * v[iii] for iii in var_rxn_inds))
+ @objective(cbm, Min, sum(rand() * v[iii] for iii in var_rxn_inds))
else
- @objective(cbmodel, Min, v[var_rxn_inds[i]])
+ @objective(cbm, Min, v[var_rxn_inds[i]])
end
- optimize!(cbmodel)
+ optimize!(cbm)
for j = 1:size(wpoints, 1)
wpoints[j, ii+1] = value(v[j])
end
@@ -56,11 +59,13 @@ function get_warmup_points(cbmodel, v, lbs, ubs; random_objective = false, numst
end
"""
- get_bound_vectors(ubconref, lbconref)
+ get_bound_vectors(opt_model)
Return Float64 vectors of the lower and upper bounds of the JuMP constraint refs.
"""
-function get_bound_vectors(lbconref, ubconref)
+function get_bound_vectors(opt_model)
+ lbconref = opt_model[:lbs]
+ ubconref = opt_model[:ubs]
lbs = zeros(length(lbconref))
for i in eachindex(lbs)
lbval = normalized_rhs(lbconref[i])
@@ -101,25 +106,24 @@ function hit_and_run(
sense = MOI.MAX_SENSE,
)
# get core optimization problem
- cbmodel, v, mb, lbcons, ubcons =
- make_optimization_model(model, optimizer, sense = sense)
+ cbm = make_optimization_model(model, optimizer, sense = sense)
+ v = cbm[:x]
if !isempty(solver_attributes) # set other attributes
for (k, v) in solver_attributes
- set_optimizer_attribute(cbmodel, k, v)
+ set_optimizer_attribute(cbm, k, v)
end
end
# set additional constraints
for (rxnid, con) in constraints
ind = model.reactions[findfirst(model.reactions, rxnid)]
- set_bound(ind, lbcons, ubcons; lb = con[1], ub = con[2])
+ set_bound(ind, cbm; lb = con[1], ub = con[2])
end
- lbs, ubs = get_bound_vectors(lbcons, ubcons) # get actual ub and lb constraints, can't use model function because the user may have changed them in the function arguments
+ lbs, ubs = get_bound_vectors(cbm) # get actual ub and lb constraints, can't use model function because the user may have changed them in the function arguments
- wpoints =
- get_warmup_points(cbmodel, v, lbcons, ubcons, random_objective = random_objective)
+ wpoints = get_warmup_points(cbm; random_objective = random_objective)
nwpts = size(wpoints, 2) # number of warmup points generated
samples = zeros(size(wpoints, 1), samplesize) # sample storage
@@ -233,12 +237,12 @@ function achr(
solver_attributes = Dict{Any,Any}(),
random_objective = false,
)
- cbmodel, v, _, ubcons, lbcons = build_cbm(model)
+ cbm, v, _, ubcons, lbcons = build_cbm(model)
- set_optimizer(cbmodel, optimizer) # choose optimizer
+ set_optimizer(cbm, optimizer) # choose optimizer
if !isempty(solver_attributes) # set other attributes
for (k, v) in solver_attributes
- set_optimizer_attribute(cbmodel, k, v)
+ set_optimizer_attribute(cbm, k, v)
end
end
@@ -250,8 +254,7 @@ function achr(
ubs, lbs = get_bound_vectors(ubcons, lbcons) # get actual ub and lb constraints
- wpoints =
- get_warmup_points(cbmodel, v, ubcons, lbcons, random_objective = random_objective)
+ wpoints = get_warmup_points(cbm, v, ubcons, lbcons, random_objective = random_objective)
nwpts = size(wpoints, 2) # number of warmup points generated
samples = zeros(size(wpoints, 1), samplesize) # sample storage
diff --git a/test/analysis/fba.jl b/test/analysis/fba.jl
index 3aecc1f119f8183b3f9e2d586668a5e33ed4051a..4e89af7b5cf72bf4813216002029669f6927571c 100644
--- a/test/analysis/fba.jl
+++ b/test/analysis/fba.jl
@@ -1,18 +1,21 @@
@testset "Flux balance analysis" begin
cp = test_simpleLP()
- (lp, x) = flux_balance_analysis(cp, GLPK.Optimizer)
+ lp = flux_balance_analysis(cp, GLPK.Optimizer)
+ x = lp[:x]
@test termination_status(lp) === MOI.OPTIMAL
sol = COBREXA.JuMP.value.(x)
@test sol ≈ [1.0, 2.0]
- (lp, x) = flux_balance_analysis(cp, Clp.Optimizer)
+ lp = flux_balance_analysis(cp, Clp.Optimizer)
+ x = lp[:x]
@test termination_status(lp) === MOI.OPTIMAL
sol = COBREXA.JuMP.value.(x)
@test sol ≈ [1.0, 2.0]
# test the maximization of the objective
cp = test_simpleLP2()
- (lp, x) = flux_balance_analysis(cp, GLPK.Optimizer)
+ lp = flux_balance_analysis(cp, GLPK.Optimizer)
+ x = lp[:x]
@test termination_status(lp) === MOI.OPTIMAL
sol = COBREXA.JuMP.value.(x)
@test sol ≈ [-1.0, 2.0]
@@ -27,7 +30,8 @@
cp = load_model(model_path, "iJR904")
expected_optimum = 0.9219480950504393
- (lp, x) = flux_balance_analysis(cp, GLPK.Optimizer)
+ lp = flux_balance_analysis(cp, GLPK.Optimizer)
+ x = lp[:x]
@test termination_status(lp) === MOI.OPTIMAL
sol = COBREXA.JuMP.value.(x)
@test objective_value(lp) ≈ expected_optimum
@@ -82,14 +86,16 @@ end
@test isapprox(producing["atp_c"]["PYK"], 1.75818, atol = 1e-3)
# set bounds
- cbm, v, mb, lbs, ubs = make_optimization_model(model, optimizer)
+ cbm = make_optimization_model(model, optimizer)
+ ubs = cbm[:ubs]
+ lbs = cbm[:lbs]
glucose_index = model[findfirst(model.reactions, "EX_glc__D_e")]
o2_index = model[findfirst(model.reactions, "EX_o2_e")]
atpm_index = model[findfirst(model.reactions, "ATPM")]
- set_bound(glucose_index, lbs, ubs; ub = -1.0, lb = -1.0)
+ set_bound(glucose_index, cbm; ub = -1.0, lb = -1.0)
@test normalized_rhs(ubs[glucose_index]) == -1.0
@test normalized_rhs(lbs[glucose_index]) == 1.0
- set_bound(o2_index, lbs, ubs; ub = 1.0, lb = 1.0)
+ set_bound(o2_index, cbm; ub = 1.0, lb = 1.0)
@test normalized_rhs(ubs[o2_index]) == 1.0
@test normalized_rhs(lbs[o2_index]) == -1.0
end
diff --git a/test/base/solver.jl b/test/base/solver.jl
index bcfce612d7af4e263510db5b421b913e9629b2de..5f7f8ab7954e6657fb8cee911e84586d2d61dd30 100644
--- a/test/base/solver.jl
+++ b/test/base/solver.jl
@@ -2,13 +2,15 @@
@testset "Solve LP" begin
cp = test_simpleLP()
optimizer = GLPK.Optimizer
- (lp, x) = optimize_model(cp, optimizer)
+ lp = optimize_model(cp, optimizer)
+ x = lp[:x]
@test termination_status(lp) === MOI.OPTIMAL
sol = COBREXA.JuMP.value.(x)
@test sol ≈ [1.0, 2.0]
optimizer = Clp.Optimizer
- (lp, x) = optimize_model(cp, optimizer)
+ lp = optimize_model(cp, optimizer)
+ x = lp[:x]
@test termination_status(lp) === MOI.OPTIMAL
sol = COBREXA.JuMP.value.(x)
@test sol ≈ [1.0, 2.0]