Commit bfb68846 authored by St. Elmo's avatar St. Elmo
Browse files

serial fva using cobramodel implemented

parent 1ff07dd5
...@@ -38,7 +38,7 @@ function fba( ...@@ -38,7 +38,7 @@ function fba(
sense = MOI.MAX_SENSE sense = MOI.MAX_SENSE
) )
# get core optimization problem # get core optimization problem
cbm, v, mb, lbcons, ubcons = makeOptimizationModel(model, optimizer) cbm, v, mb, lbcons, ubcons = makeOptimizationModel(model, optimizer, sense=sense)
# modify core optimization problem according to user specifications # modify core optimization problem according to user specifications
if !isempty(solver_attributes) # set other attributes if !isempty(solver_attributes) # set other attributes
...@@ -80,7 +80,7 @@ function fba( ...@@ -80,7 +80,7 @@ function fba(
end end
end end
@objective(cbm, Max, sum(opt_weights[i] * v[i] for i in objective_indices)) @objective(cbm, sense, sum(opt_weights[i] * v[i] for i in objective_indices))
else # use default objective else # use default objective
# automatically assigned by makeOptimizationModel # automatically assigned by makeOptimizationModel
end end
......
...@@ -96,16 +96,16 @@ fva_max, fva_min = fva(model, biomass, optimizer; solver_attributes=atts) ...@@ -96,16 +96,16 @@ fva_max, fva_min = fva(model, biomass, optimizer; solver_attributes=atts)
""" """
function fva( function fva(
model::CobraModel, model::CobraModel,
objective_rxns::Union{Reaction,Array{Reaction,1}},
optimizer; optimizer;
objective_func::Union{Reaction,Array{Reaction,1}} = Reaction[],
optimum_bound = 0.9999, optimum_bound = 0.9999,
weights = Float64[], weights = Float64[],
solver_attributes = Dict{Any,Any}(), solver_attributes = Dict{Any,Any}(),
constraints = Dict{String,Tuple{Float64,Float64}}(), constraints = Dict{String,Tuple{Float64,Float64}}(),
sense = MOI.MAX_SENSE
) )
cbm, _, _, ubcons, lbcons = build_cbm(model) # get the base constraint based model cbm, v, mb, lbcons, ubcons = makeOptimizationModel(model, optimizer, sense=sense)
set_optimizer(cbm, optimizer) # choose optimizer
if !isempty(solver_attributes) # set other attributes if !isempty(solver_attributes) # set other attributes
for (k, v) in solver_attributes for (k, v) in solver_attributes
set_optimizer_attribute(cbm, k, v) set_optimizer_attribute(cbm, k, v)
...@@ -115,35 +115,38 @@ function fva( ...@@ -115,35 +115,38 @@ function fva(
# set additional constraints # set additional constraints
for (rxnid, con) in constraints for (rxnid, con) in constraints
ind = model.reactions[findfirst(model.reactions, rxnid)] ind = model.reactions[findfirst(model.reactions, rxnid)]
set_bound(ind, ubcons, lbcons; ub = con[1], lb = con[2]) set_bound(ind, lbcons, ubcons; lb = con[1], ub = con[2])
end end
# ensure that an array of objective indices are fed in # if an objective function is supplied, modify the default objective
if typeof(objective_rxns) == Reaction if typeof(objective_func) == Reaction || !isempty(objective_func)
objective_indices = [model[objective_rxns]] # ensure that an array of objective indices are fed in
else if typeof(objective_func) == Reaction
objective_indices = [model[rxn] for rxn in objective_rxns] objective_indices = [model[objective_func]]
end else
objective_indices = [model[rxn] for rxn in objective_func]
end
if isempty(weights) if isempty(weights)
weights = ones(length(objective_indices)) weights = ones(length(objective_indices))
end end
opt_weights = zeros(length(model.reactions)) opt_weights = zeros(length(model.reactions))
# update the objective function tracker # update the objective function tracker
wcounter = 1 # don't update model objective function - silly thing to do
for i in eachindex(model.reactions) wcounter = 1
if i in objective_indices for i in eachindex(model.reactions)
model.reactions[i].objective_coefficient = weights[wcounter] if i in objective_indices
opt_weights[i] = weights[wcounter] # model.reactions[i].objective_coefficient = weights[wcounter]
wcounter += 1 opt_weights[i] = weights[wcounter]
else wcounter += 1
model.reactions[i].objective_coefficient = 0.0 # else
# model.reactions[i].objective_coefficient = 0.0
end
end end
@objective(cbm, sense, sum(opt_weights[i] * v[i] for i in objective_indices))
end end
v = all_variables(cbm)
@objective(cbm, Max, sum(opt_weights[i] * v[i] for i in objective_indices))
optimize!(cbm) optimize!(cbm)
status = ( status = (
...@@ -159,6 +162,7 @@ function fva( ...@@ -159,6 +162,7 @@ function fva(
return fva_max, fva_min return fva_max, fva_min
end end
# Now do FVA
λ = objective_value(cbm) # objective value λ = objective_value(cbm) # objective value
@constraint( @constraint(
cbm, cbm,
......
...@@ -33,14 +33,14 @@ function pfba( ...@@ -33,14 +33,14 @@ function pfba(
if typeof(optimizer) <: AbstractArray # choose optimizer if typeof(optimizer) <: AbstractArray # choose optimizer
cbm, v, mb, lbcons, ubcons = makeOptimizationModel(model, optimizer[1]) cbm, v, mb, lbcons, ubcons = makeOptimizationModel(model, optimizer[1], sense=sense)
if !isempty(solver_attributes["opt1"]) # set other attributes if !isempty(solver_attributes["opt1"]) # set other attributes
for (k, v) in solver_attributes["opt1"] for (k, v) in solver_attributes["opt1"]
set_optimizer_attribute(cbm, k, v) set_optimizer_attribute(cbm, k, v)
end end
end end
else # singe optimizer else # singe optimizer
cbm, v, mb, lbcons, ubcons = makeOptimizationModel(model, optimizer) cbm, v, mb, lbcons, ubcons = makeOptimizationModel(model, optimizer, sense=sense)
if !isempty(solver_attributes) # set other attributes if !isempty(solver_attributes) # set other attributes
for (k, v) in solver_attributes for (k, v) in solver_attributes
set_optimizer_attribute(cbm, k, v) set_optimizer_attribute(cbm, k, v)
...@@ -79,7 +79,7 @@ function pfba( ...@@ -79,7 +79,7 @@ function pfba(
end end
end end
@objective(cbm, Max, sum(opt_weights[i] * v[i] for i in objective_indices)) @objective(cbm, sense, sum(opt_weights[i] * v[i] for i in objective_indices))
else else
# objective_indices = findnz(objective(model)) # objective_indices = findnz(objective(model))
# opt_weights = ones(length(objective_indices)) # assume equal weighting, assume sense is max # opt_weights = ones(length(objective_indices)) # assume equal weighting, assume sense is max
......
# @testset "Flux variability analysis" begin @testset "Flux variability analysis" begin
# cp = test_simpleLP() cp = test_simpleLP()
# optimizer = GLPK.Optimizer optimizer = GLPK.Optimizer
# fluxes = fluxVariabilityAnalysis(cp, optimizer) fluxes = fluxVariabilityAnalysis(cp, optimizer)
# @test size(fluxes) == (2, 2) @test size(fluxes) == (2, 2)
# @test fluxes ≈ [ @test fluxes [
# 1.0 1.0 1.0 1.0
# 2.0 2.0 2.0 2.0
# ] ]
# fluxes = fluxVariabilityAnalysis(cp, [2], optimizer) fluxes = fluxVariabilityAnalysis(cp, [2], optimizer)
# @test size(fluxes) == (1, 2) @test size(fluxes) == (1, 2)
# @test fluxes == Array{Float64,2}([2 2]) @test fluxes == Array{Float64,2}([2 2])
# # a special testcase for slightly sub-optimal FVA (gamma<1) # a special testcase for slightly sub-optimal FVA (gamma<1)
# cp = LinearModel( cp = LinearModel(
# [-1.0 -1.0 -1.0], [-1.0 -1.0 -1.0],
# [0.0], [0.0],
# [1.0, 0.0, 0.0], [1.0, 0.0, 0.0],
# [0.0, 0.0, -1.0], [0.0, 0.0, -1.0],
# 1.0 * ones(3), 1.0 * ones(3),
# ["r$x" for x = 1:3], ["r$x" for x = 1:3],
# ["m1"], ["m1"],
# ) )
# fluxes = fluxVariabilityAnalysis(cp, optimizer) fluxes = fluxVariabilityAnalysis(cp, optimizer)
# @test fluxes ≈ [ @test fluxes [
# 1.0 1.0 1.0 1.0
# 0.0 0.0 0.0 0.0
# -1.0 -1.0 -1.0 -1.0
# ] ]
# fluxes = fluxVariabilityAnalysis(cp, optimizer; gamma = 0.5) fluxes = fluxVariabilityAnalysis(cp, optimizer; gamma = 0.5)
# @test fluxes ≈ [ @test fluxes [
# 0.5 1.0 0.5 1.0
# 0.0 0.5 0.0 0.5
# -1.0 -0.5 -1.0 -0.5
# ] ]
# fluxes = fluxVariabilityAnalysis(cp, optimizer; gamma = 0.0) fluxes = fluxVariabilityAnalysis(cp, optimizer; gamma = 0.0)
# @test fluxes ≈ [ @test fluxes [
# 0.0 1.0 0.0 1.0
# 0.0 1.0 0.0 1.0
# -1.0 0.0 -1.0 0.0
# ] ]
# @test isempty(fluxVariabilityAnalysis(cp, Vector{Int}(), GLPK.Optimizer)) @test isempty(fluxVariabilityAnalysis(cp, Vector{Int}(), GLPK.Optimizer))
# @test_throws DomainError fluxVariabilityAnalysis(cp, [-1], GLPK.Optimizer) @test_throws DomainError fluxVariabilityAnalysis(cp, [-1], GLPK.Optimizer)
# @test_throws DomainError fluxVariabilityAnalysis(cp, [99999999], GLPK.Optimizer) @test_throws DomainError fluxVariabilityAnalysis(cp, [99999999], GLPK.Optimizer)
# end end
# @testset "Parallel FVA" begin @testset "Parallel FVA" begin
# cp = test_simpleLP() cp = test_simpleLP()
# pids = addprocs(2, topology = :master_worker) pids = addprocs(2, topology = :master_worker)
# @everywhere using COBREXA, GLPK @everywhere using COBREXA, GLPK
# fluxes = fluxVariabilityAnalysis(cp, [1, 2], GLPK.Optimizer, pids) fluxes = fluxVariabilityAnalysis(cp, [1, 2], GLPK.Optimizer, pids)
# @test fluxes ≈ [ @test fluxes [
# 1.0 1.0 1.0 1.0
# 2.0 2.0 2.0 2.0
# ] ]
# rmprocs(pids) rmprocs(pids)
# end end
@testset "Flux variability analysis with CobraModel" begin
model = read_model(download("http://bigg.ucsd.edu/static/models/e_coli_core.json", joinpath("data", "e_coli_core.json")))
@test length(model.reactions) == 95 # read in correctly
# FVA
optimizer = Tulip.Optimizer
atts = Dict("IPM_IterationsLimit" => 400)
cons = Dict("EX_glc__D_e" => (-10.0, -10.0))
fva_max, fva_min = fva(model, biomass, optimizer, solver_attributes = atts)
fva_max2, fva_min2 =
fva(model, [biomass, pfl], optimizer, weights = [0.5, 0.5], constraints = cons)
@testset "FVA" begin
@test isapprox(fva_max["PDH"]["PDH"], 9.338922420065819, atol = 1e-6)
@test isapprox(fva_min["PDH"]["PDH"], 9.270274952732315, atol = 1e-6)
@test !isempty(fva_max2)
@test !isempty(fva_min2)
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