Unverified Commit 0e04a8bb authored by St. Elmo's avatar St. Elmo
Browse files

added constraint modification callback

parent 34b97b5c
......@@ -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
......@@ -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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment