Skip to content
Snippets Groups Projects
Unverified Commit b87478e9 authored by St. Elmo's avatar St. Elmo
Browse files

small updates

parent 0e04a8bb
No related branches found
No related tags found
No related merge requests found
...@@ -89,6 +89,7 @@ sense = MOI.MAX_SENSE ...@@ -89,6 +89,7 @@ sense = MOI.MAX_SENSE
end end
end end
<<<<<<< Updated upstream
# if an objective function is supplied, modify the default objective # if an objective function is supplied, modify the default objective
if typeof(objective_func) == Reaction || !isempty(objective_func) if typeof(objective_func) == Reaction || !isempty(objective_func)
# ensure that an array of objective indices are fed in # ensure that an array of objective indices are fed in
...@@ -123,10 +124,32 @@ sense = MOI.MAX_SENSE ...@@ -123,10 +124,32 @@ sense = MOI.MAX_SENSE
optimize!(cbm) optimize!(cbm)
=======
function flux_balance_analysis_vec(args)
cbm = flux_balance_analysis(args...)
status = ( status = (
termination_status(cbm) == MOI.OPTIMAL || termination_status(cbm) == MOI.OPTIMAL ||
termination_status(cbm) == MOI.LOCALLY_SOLVED termination_status(cbm) == MOI.LOCALLY_SOLVED
) )
if status
return value.(cbm[:x])
else
@warn "Optimization issues occurred."
return Vector{Float64}[]
end
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 if status
return map_fluxes(v, model) return map_fluxes(v, model)
...@@ -134,4 +157,13 @@ sense = MOI.MAX_SENSE ...@@ -134,4 +157,13 @@ sense = MOI.MAX_SENSE
@warn "Optimization issues occurred." @warn "Optimization issues occurred."
return Dict{String,Float64}() return Dict{String,Float64}()
end end
=======
if status
return Dict(zip(reactions(model), value.(cbm[:x])))
else
@warn "Optimization issues occurred."
return Dict{String, Float64}()
end
>>>>>>> Stashed changes
end end
...@@ -38,7 +38,7 @@ function flux_variability_analysis( ...@@ -38,7 +38,7 @@ function flux_variability_analysis(
end end
optimization_model = flux_balance_analysis(model, optimizer) optimization_model = flux_balance_analysis(model, optimizer)
Z0 = JuMP.objective_value(optimization_model) Z0 = COBREXA.JuMP.objective_value(optimization_model)
optimization_model = nothing # we won't need this one anymore, so free the memory optimization_model = nothing # we won't need this one anymore, so free the memory
# store a JuMP optimization model at all workers # store a JuMP optimization model at all workers
...@@ -92,7 +92,7 @@ Internal helper function for adding constraints to a model. Exists mainly ...@@ -92,7 +92,7 @@ Internal helper function for adding constraints to a model. Exists mainly
because for avoiding namespace problems on remote workers. because for avoiding namespace problems on remote workers.
""" """
function _FVA_add_constraint(model, c, x, Z0, gamma) function _FVA_add_constraint(model, c, x, Z0, gamma)
JuMP.@constraint(model, c' * x gamma * Z0) COBREXA.JuMP.@constraint(model, c' * x gamma * Z0)
end end
""" """
...@@ -103,11 +103,11 @@ namespace problems. ...@@ -103,11 +103,11 @@ namespace problems.
""" """
function _FVA_optimize_reaction(model, rid) function _FVA_optimize_reaction(model, rid)
sense = rid > 0 ? MOI.MAX_SENSE : MOI.MIN_SENSE sense = rid > 0 ? MOI.MAX_SENSE : MOI.MIN_SENSE
var = JuMP.all_variables(model)[abs(rid)] var = COBREXA.JuMP.all_variables(model)[abs(rid)]
JuMP.@objective(model, sense, var) COBREXA.JuMP.@objective(model, sense, var)
JuMP.optimize!(model) COBREXA.JuMP.optimize!(model)
return JuMP.objective_value(model) return COBREXA.JuMP.objective_value(model)
end end
""" """
...@@ -134,20 +134,60 @@ fva_max, fva_min = fva(model, biomass, optimizer; solver_attributes=atts) ...@@ -134,20 +134,60 @@ fva_max, fva_min = fva(model, biomass, optimizer; solver_attributes=atts)
function fva( function fva(
model::CobraModel, model::CobraModel,
optimizer; optimizer;
objective_func::Union{Reaction,Array{Reaction,1}} = Reaction[],
optimum_bound = 0.9999, optimum_bound = 0.9999,
modifications = [(model, opt_model) -> nothing] weights = Float64[],
solver_attributes = Dict{Any,Any}(),
constraints = Dict{String,Tuple{Float64,Float64}}(),
sense = MOI.MAX_SENSE,
) )
<<<<<<< Updated upstream
cbm = make_optimization_model(model, optimizer, sense = sense)
v = cbm[:x]
=======
# get core optimization problem # get core optimization problem
cbm = make_optimization_model(model, optimizer) cbm = make_optimization_model(model, optimizer)
v = cbm[:x] >>>>>>> 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
# apply callbacks - user can also just put in a function # if an objective function is supplied, modify the default objective
if typeof(modifications) <: Vector if typeof(objective_func) == Reaction || !isempty(objective_func)
for mod in modifications # ensure that an array of objective indices are fed in
mod(model, cbm) if typeof(objective_func) == Reaction
end objective_indices = [model[objective_func]]
else else
modifications(model, cbm) 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 end
optimize!(cbm) optimize!(cbm)
...@@ -166,6 +206,8 @@ function fva( ...@@ -166,6 +206,8 @@ function fva(
end end
# Now do FVA # Now do FVA
v = cbm[:x]
λ = objective_value(cbm) # objective value λ = objective_value(cbm) # objective value
@constraint( @constraint(
cbm, cbm,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment