From 0e04a8bb7dccf252df6377edb0e6e3af6c9364bc Mon Sep 17 00:00:00 2001
From: stelmo <stelmozors@gmail.com>
Date: Thu, 1 Apr 2021 14:47:49 +0200
Subject: [PATCH] added constraint modification callback
---
src/analysis/fba.jl | 95 ++++++++++++++++++++++++++++++-------------
src/base/utilities.jl | 60 ---------------------------
2 files changed, 66 insertions(+), 89 deletions(-)
diff --git a/src/analysis/fba.jl b/src/analysis/fba.jl
index ad8a5c748..7b9beed1e 100644
--- a/src/analysis/fba.jl
+++ b/src/analysis/fba.jl
@@ -22,9 +22,10 @@ Arguments are passed to [`flux_balance_analysis`](@ref).
"""
function flux_balance_analysis_vec(args...)::Union{Vector{Float64},Nothing}
optmodel = flux_balance_analysis(args...)
-
+ vars = optmodel[:x]
+
termination_status(optmodel) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
- value.(optmodel[:x])
+ value.(vars)
end
"""
@@ -43,7 +44,7 @@ function flux_balance_analysis_dict(
end
"""
- flux_balance_analysis(model::CobraModel, optimizer; objective_func::Union{Reaction, Array{Reaction, 1}}=Reaction[], weights=Float64[], solver_attributes=Dict{Any, Any}(), constraints=Dict{String, Tuple{Float64,Float64}}())
+ fba(model::CobraModel, optimizer; objective_func::Union{Reaction, Array{Reaction, 1}}=Reaction[], weights=Float64[], solver_attributes=Dict{Any, Any}(), constraints=Dict{String, Tuple{Float64,Float64}}())
Run flux balance analysis (FBA) on the `model` optionally specifying `objective_rxn(s)` and their `weights` (empty `weights` mean equal weighting per reaction).
Optionally also specify any additional flux constraints with `constraints`, a dictionary mapping reaction `id`s to tuples of (lb, ub) flux constraints.
@@ -61,40 +62,76 @@ biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
sol = fba(model, biomass, optimizer; solver_attributes=atts)
```
"""
-function flux_balance_analysis(
+function fba(
model::CobraModel,
optimizer;
- modifications::Union{Vector, Function} = [(model, opt_model)->nothing]
+ modifications = [(model, opt_model)->nothing]
)
+
+objective_func::Union{Reaction,Array{Reaction,1}} = Reaction[]
+weights = Float64[]
+solver_attributes = Dict{Any,Any}()
+sense = MOI.MAX_SENSE
+
# get core optimization problem
- cbm = make_optimization_model(model, optimizer)
-
- # apply callbacks - user can also just put in a function
- if typeof(modifications) <: Vector
- for mod in modifications
- mod(model, cbm)
- end
- else
- modifications(model, cbm)
+ cbm = make_optimization_model(model, optimizer, sense = sense)
+ v = cbm[:x] # fluxes
+
+ # apply callbacks
+ for mod in modifications
+ mod(model, cbm)
end
- optimize!(cbm)
-
- return cbm
-end
+ # modify core optimization problem according to user specifications
+ if !isempty(solver_attributes) # set other attributes
+ for (k, val) in solver_attributes
+ set_optimizer_attribute(cbm, k, val)
+ end
+ end
-function flux_balance_analysis_vec()
+ # if an objective function is supplied, modify the default objective
+ if typeof(objective_func) == Reaction || !isempty(objective_func)
+ # ensure that an array of objective indices are fed in
+ if typeof(objective_func) == Reaction
+ objective_indices = [model[objective_func]]
+ else
+ objective_indices = [model[rxn] for rxn in objective_func]
+ end
-end
+ if isempty(weights)
+ weights = ones(length(objective_indices))
+ end
+ opt_weights = zeros(length(model.reactions))
-status = (
- termination_status(cbm) == MOI.OPTIMAL ||
- termination_status(cbm) == MOI.LOCALLY_SOLVED
-)
+ # update the objective function tracker
+ # don't update model objective function - silly thing to do
+ wcounter = 1
+ for i in eachindex(model.reactions)
+ if i in objective_indices
+ # model.reactions[i].objective_coefficient = weights[wcounter]
+ opt_weights[i] = weights[wcounter]
+ wcounter += 1
+ # else
+ # model.reactions[i].objective_coefficient = 0.0
+ end
+ end
+
+ @objective(cbm, sense, sum(opt_weights[i] * v[i] for i in objective_indices))
+ else # use default objective
+ # automatically assigned by make_optimization_model
+ end
+
+ optimize!(cbm)
+
+ status = (
+ termination_status(cbm) == MOI.OPTIMAL ||
+ termination_status(cbm) == MOI.LOCALLY_SOLVED
+ )
-if status
- return map_fluxes(cbm[:v], model)
-else
- @warn "Optimization issues occurred."
- return Dict{String,Float64}()
+ if status
+ return map_fluxes(v, model)
+ else
+ @warn "Optimization issues occurred."
+ return Dict{String,Float64}()
+ end
end
diff --git a/src/base/utilities.jl b/src/base/utilities.jl
index 0ae3f14c9..e2d6d7233 100644
--- a/src/base/utilities.jl
+++ b/src/base/utilities.jl
@@ -46,63 +46,3 @@ function modify_constraint(reaction::Reaction, lb, ub)
set_bound(ind, opt_model, lb=lb, ub=ub)
end
end
-
-"""
- modify_sense(sense)
-
- Modify the objective sense.
-Allowed options are MOI.MAX_SENSE and MOI.MIN_SENSE.
-"""
-function modify_sense(sense)
- (model, opt_model) -> begin
- set_objective_sense(opt_model, sense)
- end
-end
-
-"""
-modify_solver_attributes(option_key, option_value)
-
-"""
-function modify_solver_attributes(option_key, option_value)
- (model, opt_model) -> begin
- set_optimizer_attribute(opt_model, option_key, option_value)
- end
-end
-
-"""
- modify_objective(objective_function::Union{Reaction,Array{Reaction,1}}; objective_weights=[], sense=MOI.MAX_SENSE)
-
-Modify the objective function and optionally specify weights
-(useful if more than one reaction is set in `objective_rxns`)
-and the sense of the objective.
-
-# Example
-"""
-function modify_objective(objective_function::Union{Reaction,Array{Reaction,1}}; objective_weights=[], sense=MOI.MAX_SENSE)
- (model, opt_model) -> begin
- v = opt_model[:x]
-
- # ensure that an array of objective indices are fed in
- if typeof(objective_function) == Reaction
- objective_indices = [model[objective_function]]
- else
- objective_indices = [model[rxn] for rxn in objective_function]
- end
-
- if isempty(objective_weights)
- weights = ones(length(objective_indices))
- end
- opt_weights = zeros(length(model.reactions))
-
- # update the objective function tracker
- wcounter = 1
- for i in eachindex(model.reactions)
- if i in objective_indices
- opt_weights[i] = weights[wcounter]
- wcounter += 1
- end
- end
-
- @objective(cbm, MOI.MAX_SENSE, sum(opt_weights[i] * v[i] for i in objective_indices))
- end
-end
--
GitLab