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

fixed the chaos :)

parent b87478e9
......@@ -22,10 +22,9 @@ 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.(vars)
JuMP.termination_status(optmodel) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
value.(optmodel[:x])
end
"""
......@@ -44,126 +43,80 @@ function flux_balance_analysis_dict(
end
"""
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}}())
flux_balance_analysis(model::StandardModel, optimizer; modifications)
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.
Run flux balance analysis (FBA) on the `model` optionally specifying `modifications` to the problem.
These modifications can be entered as an array of modifications, or a single modification.
Leave this keyword argument out if you do not want to modify the problem.
See [`modify_constraint`](@ref), [`modify_solver_attribute`](@ref),[`modify_objective`](@ref), and [`modify_sense`](@ref)
for possible modifications.
Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.
The `solver_attributes` can also be specified in the form of a dictionary where each (key, value) pair will be passed to `set_optimizer_attribute(cbmodel, key, value)`.
This function builds the optimization problem from the model, and hence uses the constraints implied by the model object.
Returns a dictionary of reaction `id`s mapped to fluxes if solved successfully, otherwise an empty dictionary.
Returns a solved JuMP model.
# Example
```
optimizer = Gurobi.Optimizer
atts = Dict("OutputFlag" => 0)
model = CobraTools.read_model("iJO1366.json")
biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
sol = fba(model, biomass, optimizer; solver_attributes=atts)
model = CobraTools.read_model("e_coli_core.json")
biomass = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")
solved_model = fba(model, optimizer; modifications=[modify_objective(biomass)])
```
"""
function fba(
model::CobraModel,
function flux_balance_analysis(
model::StandardModel,
optimizer;
modifications = [(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, sense = sense)
v = cbm[:x] # fluxes
cbm = make_optimization_model(model, optimizer)
# apply callbacks
for mod in modifications
mod(model, 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)
if typeof(modifications) <: Vector # many modifications
for mod in modifications
mod(model, cbm)
end
else # single modification
modifications(model, cbm)
end
<<<<<<< Updated upstream
# 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
JuMP.optimize!(cbm)
if isempty(weights)
weights = ones(length(objective_indices))
end
opt_weights = zeros(length(model.reactions))
# 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
return cbm
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
"""
flux_balance_analysis_vec(model::StandardModel, optimizer; modifications)
optimize!(cbm)
Perform flux balance analysis on `model` using `optimizer`.
Returns a vector fluxes in the order supplied in `model`.
Calls [`flux_balance_analysis`](@ref) internally.
"""
function flux_balance_analysis_vec(
model::StandardModel,
optimizer;
modifications = [(model, opt_model) -> nothing],
)
cbm = flux_balance_analysis(model, optimizer; modifications = modifications)
=======
function flux_balance_analysis_vec(args)
cbm = flux_balance_analysis(args...)
status = (
termination_status(cbm) == MOI.OPTIMAL ||
termination_status(cbm) == MOI.LOCALLY_SOLVED
)
if status
return value.(cbm[:x])
else
@warn "Optimization issues occurred."
return Vector{Float64}[]
end
JuMP.termination_status(cbm) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
return value.(cbm[:x])
end
function flux_balance_analysis_dict(model, args)
cbm = flux_balance_analysis(model, args...)
>>>>>>> Stashed changes
status = (
termination_status(cbm) == MOI.OPTIMAL ||
termination_status(cbm) == MOI.LOCALLY_SOLVED
)
<<<<<<< Updated upstream
if status
return map_fluxes(v, model)
else
@warn "Optimization issues occurred."
return Dict{String,Float64}()
end
=======
"""
flux_balance_analysis_dict(model::StandardModel, optimizer; modifications)
Perform flux balance analysis on `model` using `optimizer`.
Returns a dictionary mapping reaction `id`s to fluxes.
Calls [`flux_balance_analysis`](@ref) internally.
"""
function flux_balance_analysis_dict(
model::StandardModel,
optimizer;
modifications = [(model, opt_model) -> nothing],
)
cbm = flux_balance_analysis(model, optimizer; modifications = modifications)
JuMP.termination_status(cbm) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
if status
return Dict(zip(reactions(model), value.(cbm[:x])))
else
@warn "Optimization issues occurred."
return Dict{String, Float64}()
end
>>>>>>> Stashed changes
return Dict(zip(reactions(model), value.(cbm[:x])))
end
......@@ -45,8 +45,13 @@ function flux_variability_analysis(
save_model = :(
begin
optmodel = COBREXA.make_optimization_model($model, $optimizer)
x = optmodel[:x]
COBREXA._FVA_add_constraint(optmodel, $(objective(model)), x, $Z0, $gamma)
COBREXA._FVA_add_constraint(
optmodel,
$(objective(model)),
optmodel[:x],
$Z0,
$gamma,
)
optmodel
end
)
......@@ -111,120 +116,60 @@ function _FVA_optimize_reaction(model, rid)
end
"""
fva(model::CobraModel, objective_rxns::Union{Reaction, Array{Reaction, 1}}, optimizer; optimum_bound=0.9999, weights=Float64[], solver_attributes=Dict{Any, Any}(), constraints=Dict{String, Tuple{Float64,Float64}}())
fva(model::StandardModel, optimizer; optimum_bound=0.9999, modifications)
Run flux variability analysis (FVA) on the `model` with `objective_rxn(s)` and optionally specifying their `weights` (empty `weights` mean equal weighting per reaction).
It runs fba on the model once to determine the optimum of the objective.
Optionally also specify any additional flux constraints with `constraints`, a dictionary mapping reaction `id`s to tuples of (ub, lb) flux constraints.
The model is then constrained to produce objective flux bounded by `optimum_bound` from below (set to slightly less than 1.0 for stability) and each flux in the model sequentially minimized and maximized.
Run flux variability analysis (FVA) on the `model` (of type `StandardModel`).
Optionally specifying problem modifications like in [`flux_balance_analysis`](@ref).
This algorithm runs FBA on the model to determine the optimum of the objective.
This optimum then constrains subsequent problems, where `optimum_bound` can be used to relax this constraint as a fraction of the FBA optimum.
Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.
The `solver_attributes` can also be specified in the form of a dictionary where each (key, value) pair will be passed to `set_optimizer_attribute(cbmodel, key, value)`.
This function builds the optimization problem from the model, and hence uses the constraints implied by the model object.
Returns two dictionaries (`fva_max` and `fva_min`) that each reaction `id`s to dictionaries of the resultant flux distributions (if solved successfully) when that `id` is optimized.
Note, this function only runs serially.
Consider using a different model type for parallel implementations.
Returns two dictionaries (`fva_max` and `fva_min`) that map each reaction `id` to dictionaries of the resultant flux distributions when that `id` is optimized.
# Example
```
optimizer = Gurobi.Optimizer
atts = Dict("OutputFlag" => 0)
model = Cobraread_model("iJO1366.json")
biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
fva_max, fva_min = fva(model, biomass, optimizer; solver_attributes=atts)
```
"""
function fva(
model::CobraModel,
function flux_variability_analysis(
model::StandardModel,
optimizer;
objective_func::Union{Reaction,Array{Reaction,1}} = Reaction[],
optimum_bound = 0.9999,
weights = Float64[],
solver_attributes = Dict{Any,Any}(),
constraints = Dict{String,Tuple{Float64,Float64}}(),
sense = MOI.MAX_SENSE,
modifications = [(model, opt_model) -> nothing],
)
<<<<<<< Updated upstream
cbm = make_optimization_model(model, optimizer, sense = sense)
v = cbm[:x]
=======
# get core optimization problem
cbm = make_optimization_model(model, optimizer)
>>>>>>> Stashed changes
if !isempty(solver_attributes) # set other attributes
for (k, v) in solver_attributes
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, cbm; lb = con[1], ub = con[2])
end
# 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
if isempty(weights)
weights = ones(length(objective_indices))
end
opt_weights = zeros(length(model.reactions))
# 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))
end
optimize!(cbm)
status = (
termination_status(cbm) == MOI.OPTIMAL ||
termination_status(cbm) == MOI.LOCALLY_SOLVED
)
# Run FBA
cbm = flux_balance_analysis(model, optimizer; modifications=modifications)
fva_min = Dict{String,Dict{String,Float64}}()
fva_max = Dict{String,Dict{String,Float64}}()
JuMP.termination_status(cbm) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || (return nothing, nothing)
if !status
@warn "Error getting the initial optimum, aborting "
return fva_max, fva_min
end
fva_min = Dict{String,Union{Nothing, Dict{String,Float64}}}()
fva_max = Dict{String,Union{Nothing, Dict{String,Float64}}}()
# Now do FVA
v = cbm[:x]
λ = objective_value(cbm) # objective value
@constraint(
λ = JuMP.objective_value(cbm) # objective value
JuMP.@constraint(
cbm,
optimum_bound * λ <= sum(opt_weights[i] * v[i] for i in objective_indices) <= λ
optimum_bound * λ <= JuMP.objective_function(cbm) <= λ * 1.0/optimum_bound # in case there is a negative bound and the problem has min sense
) # for stability
for i = 1:length(v)
@objective(cbm, Max, v[i])
optimize!(cbm)
JuMP.@objective(cbm, Max, v[i])
JuMP.optimize!(cbm)
status = (
termination_status(cbm) == MOI.OPTIMAL ||
termination_status(cbm) == MOI.LOCALLY_SOLVED
JuMP.termination_status(cbm) == MOI.OPTIMAL ||
JuMP.termination_status(cbm) == MOI.LOCALLY_SOLVED
)
if status
fva_max[model.reactions[i].id] = map_fluxes(v, model)
fva_max[model.reactions[i].id] = Dict(zip(reactions(model), value.(cbm[:x])))
else
@warn "Error maximizing index: $i with error $(termination_status(cbm))"
fva_max[model.reactions[i].id] = nothing
end
@objective(cbm, Min, v[i])
......@@ -234,9 +179,10 @@ function fva(
termination_status(cbm) == MOI.LOCALLY_SOLVED
)
if status
fva_min[model.reactions[i].id] = map_fluxes(v, model)
fva_min[model.reactions[i].id] = Dict(zip(reactions(model), value.(cbm[:x])))
else
@warn "Error minimizing index: $i with error $(termination_status(cbm))"
fva_min[model.reactions[i].id] = nothing
end
end
......
"""
pfba(model::CobraModel, optimizer; objective_func::Union{Reaction, Array{Reaction, 1}}=Reaction[], weights=Float64[], solver_attributes=Dict{Any, Any}(), constraints=Dict{String, Tuple{Float64,Float64}}(), sense=MOI.MAX_SENSE)
pfba(model::StandardModel, optimizer; modifications, qp_solver, qp_solver_attributes)
Run parsimonious flux balance analysis (pFBA) on the `model` with `objective_rxn(s)` and optionally specifying their `weights` (empty `weights` mean equal weighting per reaction) for the initial FBA problem.
Run parsimonious flux balance analysis (pFBA) on the `model`.
Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.
Optionally also specify any additional flux constraints with `constraints`, a dictionary mapping reaction `id`s to tuples of (lb, ub) flux constraints.
When `optimizer` is an array of optimizers, e.g. `[opt1, opt2]`, then `opt1` is used to solve the FBA problem, and `opt2` is used to solve the QP problem.
This strategy is useful when the QP solver is not good at solving the LP problem.
The `solver_attributes` can also be specified in the form of a dictionary where each (key, value) pair will be passed to `set_optimizer_attribute(cbmodel, k, v)`.
If more than one solver is specified in `optimizer`, then `solver_attributes` must be a dictionary of dictionaries with keys "opt1" and "opt2", e.g. Dict("opt1" => Dict{Any, Any}(),"opt2" => Dict{Any, Any}()).
This function builds the optimization problem from the model, and hence uses the constraints implied by the model object.
Returns a dictionary of reaction `id`s mapped to fluxes if solved successfully, otherwise an empty dictionary.
Optionally, specify problem `modifications` as in [`flux_balance_analysis`](@ref).
Also, `qp_solver` can be set to be different from `optimizer`, where the latter is then the LP optimizer only.
Also note that `qp_solver_attributes` is meant to set the attributes for the `qp_solver` and NOT `modifications`.
Any solver attributes changed in `modifications` will only affect he LP solver.
This function automatically relaxes the constraint that the FBA solution objective matches the pFBA solution.
This is iteratively relaxed like 1.0, 0.999999, 0.9999, etc. of the bound until 0.99 when the function fails and returns nothing.
Return a solved pFBA model.
# Example
```
......@@ -20,116 +20,67 @@ biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
sol = pfba(model, biomass, optimizer; solver_attributes=atts)
```
"""
function pfba(
model::CobraModel,
optimizer;
modifications = [(model, opt_model) -> nothing]
)
## FBA ################################################
function parsimonious_flux_balance_analysis(model::StandardModel, optimizer; modifications = [(model, opt_model) -> nothing], qp_solver = [(model, opt_model) -> nothing], qp_solver_attributes=[(model, opt_model) -> nothing])
# Run FBA
cbm = flux_balance_analysis(model, optimizer)#; modifications=modifications)
JuMP.termination_status(cbm) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing # FBA failed
# Run pFBA
λ = objective_value(cbm)
old_objective = JuMP.objective_function(cbm)
if typeof(optimizer) <: AbstractArray # choose optimizer
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 = 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)
end
qp_solver(model, cbm) # change the solver if specified, otherwise default argument does nothing
if typeof(qp_solver_attributes) <: Vector # many modifications
for mod in qp_solver_attributes
mod(model, cbm)
end
else # single modification
qp_solver_attributes(model, cbm)
end
# set additional constraints
for (rxnid, con) in constraints
ind = model.reactions[findfirst(model.reactions, rxnid)]
set_bound(ind, cbm; lb = con[1], ub = con[2])
end
# check if default objective should be used
if typeof(objective_func) == Reaction || !isempty(objective_func)
# check if 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
if isempty(weights)
weights = ones(length(objective_indices))
end
opt_weights = zeros(length(model.reactions))
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
v = cbm[:x] # fluxes
@objective(cbm, Min, sum(dot(v, v)))
@objective(cbm, sense, sum(opt_weights[i] * v[i] for i in objective_indices))
else
# objective_indices = findnz(objective(model))
# opt_weights = ones(length(objective_indices)) # assume equal weighting, assume sense is max
for lbconval in [1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99] # sequentially relax bound for stability
@constraint(
cbm,
pfbacon,
lbconval * λ <= old_objective <= λ * 1.0/lbconval # in case of negative objective and minimized fba
)
optimize!(cbm)
JuMP.termination_status(cbm) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] && break
COBREXA.JuMP.delete(cbm, pfbacon)
end
optimize!(cbm)
fba_status = (
termination_status(cbm) == MOI.OPTIMAL ||
termination_status(cbm) == MOI.LOCALLY_SOLVED
)
## pFBA ###############################################
λ = objective_value(cbm)
if typeof(optimizer) <: AbstractArray # choose optimizer
set_optimizer(cbm, optimizer[2])
if !isempty(solver_attributes["opt2"]) # set other attributes
for (k, v) in solver_attributes["opt2"]
set_optimizer_attribute(cbm, k, v)
end
end
end
JuMP.termination_status(cbm) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing # pFBA failed
@constraint(
cbm,
pfbacon,
λ <= sum(opt_weights[i] * v[i] for i in objective_indices) <= λ
)
@objective(cbm, Min, sum(dot(v, v)))
return cbm
end
optimize!(cbm)
"""
parsimonious_flux_balance_analysis_dict(model::StandardModel, optimizer; modifications, qp_solver, qp_solver_attributes)
for lbconval in [0.999999, 0.99999, 0.9999, 0.999, 0.99] # relax bound for stability
if termination_status(cbm) == MOI.OPTIMAL ||
termination_status(cbm) == MOI.LOCALLY_SOLVED # try to relax bound if failed optimization
break
else
COBREXA.JuMP.delete(cbm, pfbacon)
@constraint(cbm, lbconval * λ <= sum(v[i] for i in objective_indices) <= λ)
optimize!(cbm)
end
end
Perform parsimonious flux balance analysis on `model` using `optimizer`.
Returns a vector of fluxes in the same order as the reactions in `model`.
Calls [`parsimonious_flux_balance_analysis_dict`](@ref) internally.
"""
function parsimonious_flux_balance_analysis_vec(model::StandardModel, optimizer; modifications = [(model, opt_model) -> nothing], qp_solver = [(model, opt_model) -> nothing], qp_solver_attributes=[(model, opt_model) -> nothing])
cbm = parsimonious_flux_balance_analysis(model, optimizer; modifications=modifications, qp_solver=qp_solver, qp_solver_attributes=qp_solver_attributes)
JuMP.termination_status(cbm) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
return value.(cbm[:x])
end
pfba_status = (
termination_status(cbm) == MOI.OPTIMAL ||
termination_status(cbm) == MOI.LOCALLY_SOLVED
)
"""
parsimonious_flux_balance_analysis_dict(model::StandardModel, optimizer; modifications, qp_solver, qp_solver_attributes)
if fba_status && pfba_status
return map_fluxes(v, model)
else
@warn "Optimization issues occurred."
return Dict{String,Float64}() # something went wrong
end
Perform parsimonious flux balance analysis on `model` using `optimizer`.
Returns a dictionary mapping reaction `id`s to fluxes.
Calls [`parsimonious_flux_balance_analysis_dict`](@ref) internally.
"""
function parsimonious_flux_balance_analysis_dict(model::StandardModel, optimizer; modifications = [(model, opt_model) -> nothing], qp_solver = [(model, opt_model) -> nothing], qp_solver_attributes=[(model, opt_model) -> nothing])
cbm = parsimonious_flux_balance_analysis(model, optimizer; modifications=modifications, qp_solver=qp_solver, qp_solver_attributes=qp_solver_attributes)
JuMP.termination_status(cbm) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
return Dict(zip(reactions(model), value.(cbm[:x])))
end
"""
map_fluxes(v, model::CobraModel)
Map fluxes from an optimization problem (`v`) to rxns in a model.
`v` can be a JuMP object (fluxes) or an array of Float64 fluxes.
Assumes they are in order of `model.reactions`, which they should be since the optimization problem is constructed from the model.
"""
function map_fluxes(v::Array{Float64,1}, model::CobraModel)
rxndict = Dict{String,Float64}()
for i in eachindex(model.reactions)
rxndict[model.reactions[i].id] = v[i]
end
return rxndict
end
function map_fluxes(v::Array{VariableRef,1}, model::CobraModel)
rxndict = Dict{String,Float64}()
for i in eachindex(model.reactions)
rxndict[model.reactions[i].id] = value(v[i])
end
return rxndict
end
"""
atom_exchange(flux_dict::Dict{String, Float64}, model::CobraModel)
atom_exchange(flux_dict::Dict{String, Float64}, model::StandardModel)
Return a dictionary mapping the flux of atoms across the boundary of the model given `flux_dict` of reactions in `model`.
Here `flux_dict` is a mapping of reaction `id`s to fluxes, e.g. from FBA.
"""
function atom_exchange(flux_dict::Dict{String,Float64}, model::CobraModel)
function atom_exchange(flux_dict::Dict{String,Float64}, model::StandardModel)
atom_flux = Dict{String,Float64}()
for (rxnid, flux) in flux_dict
if startswith(rxnid, "EX_") || startswith(rxnid, "DM_") # exchange, demand reaction
......@@ -101,11 +78,11 @@ function exchange_reactions(
end
"""
metabolite_fluxes(fluxdict::Dict{String, Float64}, model::CobraModel)
metabolite_fluxes(fluxdict::Dict{String, Float64}, model::StandardModel)