diff --git a/src/analysis/dfba.jl b/src/analysis/dfba.jl
deleted file mode 100644
index 464090415c47109523e91779d4f40e19495c9cf1..0000000000000000000000000000000000000000
--- a/src/analysis/dfba.jl
+++ /dev/null
@@ -1 +0,0 @@
-# TODO
diff --git a/src/analysis/fba.jl b/src/analysis/fba.jl
index 3e25aa0fe6f5328f7e3df35f5e15dd65eb72e9ea..bc03a097f96ec8aff9c27f4d11c45df04fa0b975 100644
--- a/src/analysis/fba.jl
+++ b/src/analysis/fba.jl
@@ -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
diff --git a/src/analysis/fva.jl b/src/analysis/fva.jl
index 17794201f378dc2b9df0d8913531fe592f333282..19d4db4cb3a78d93747960b3588460c4e06ef256 100644
--- a/src/analysis/fva.jl
+++ b/src/analysis/fva.jl
@@ -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
diff --git a/src/analysis/pfba.jl b/src/analysis/pfba.jl
index e23f423115ae546e7ce71c40179d7984f15e2aa2..5b611222326055e371fbf03fa13e098aa0aa3fed 100644
--- a/src/analysis/pfba.jl
+++ b/src/analysis/pfba.jl
@@ -1,15 +1,15 @@
"""
- 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
diff --git a/src/analysis/fluxes.jl b/src/analysis/utilities.jl
similarity index 50%
rename from src/analysis/fluxes.jl
rename to src/analysis/utilities.jl
index 4d9db2f15bced0de35c106955803722b80530983..c23b9988c33c0d65a4b6db8184061cffe1aad340 100644
--- a/src/analysis/fluxes.jl
+++ b/src/analysis/utilities.jl
@@ -1,33 +1,10 @@
"""
- 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)
Return two dictionaries of metabolite `id`s mapped to reactions that consume or produce them given the flux distribution supplied in `fluxdict`.
"""
-function metabolite_fluxes(fluxdict::Dict{String,Float64}, model::CobraModel)
+function metabolite_fluxes(fluxdict::Dict{String,Float64}, model::StandardModel)
S = Array(stoichiometry(model))
met_flux = Dict{String,Float64}()
rxnids = reactions(model)
@@ -134,3 +111,110 @@ function metabolite_fluxes(fluxdict::Dict{String,Float64}, model::CobraModel)
end
return consuming, producing
end
+
+"""
+ set_bound(index, optimization_model; ub=1000, lb=-1000)
+Helper function to set the bounds of variables.
+The JuMP `set_normalized_rhs` function is a little confusing, so this function simplifies setting
+constraints. Just supply the constraint `index` and the model and they will be set to `ub` and `lb`.
+"""
+function set_bound(vind, opt_model; ub = 1000, lb = -1000)
+ if lb <= 0
+ set_normalized_rhs(opt_model[:lbs][vind], abs(lb))
+ else
+ set_normalized_rhs(opt_model[:lbs][vind], -abs(lb))
+ end
+ set_normalized_rhs(opt_model[:ubs][vind], ub)
+end
+
+"""
+ modify_constraint(reaction::Reaction, lb, ub)
+
+Modify constraints of model reaction.
+"""
+function modify_constraint(reaction::Reaction, lb, ub)
+ (model, opt_model) -> begin
+ ind = model.reactions[reaction]
+ set_bound(ind, opt_model, lb = lb, ub = ub)
+ end
+end
+
+"""
+ modify_solver_attribute(option_key, option_val)
+
+Modify a solver attribute. These attributes are solver specific,
+refer the either JuMP or the solver you are using's documentation.
+"""
+function modify_solver_attribute(option_key, option_val)
+ (model, opt_model) -> begin
+ JuMP.set_optimizer_attribute(opt_model, option_key, option_val)
+ end
+end
+
+"""
+ modify_sense(objective_sense)
+
+Modify the objective sense.
+Possible arguments are `MOI.MAX_SENSE` and `MOI.MIN_SENSE`.
+Note, [`modify_objective`](@ref) sets the sense of the objective,
+so it doesn't make sense to use this function AND [`modify_objective`](@ref) simultaneously.
+"""
+function modify_sense(objective_sense)
+ (model, opt_model) -> begin
+ JuMP.set_objective_sense(opt_model, objective_sense)
+ end
+end
+
+"""
+ modify_objective(objective_functions::Union{Reaction, Array{Reaction, 1}}; weights=[], sense=MOI.MAX_SENSE)
+
+Callback function to modify the objective function used in a constraint based analysis function.
+`objective_functions` can be a single reaction or an array of reactions (of type `Reaction`).
+Optionally specify their `weights` and the sense of the new objective (`MOI.MAX_SENSE`, `MOI.MIN_SENSE`).
+Note, this function sets the sense of the objective.
+"""
+function modify_objective(
+ objective_functions::Union{Reaction,Array{Reaction,1}};
+ weights = [],
+ sense = MOI.MAX_SENSE,
+)
+ (model, opt_model) -> begin
+
+ # Construct objective_indices array
+ if typeof(objective_functions) == Reaction
+ objective_indices = [model[objective_functions]]
+ else
+ objective_indices = [model[rxn] for rxn in objective_functions]
+ end
+
+ # Initialize weights
+ opt_weights = zeros(length(model.reactions))
+
+ isempty(weights) && (weights = ones(length(objective_indices))) # equal weights
+
+ wcounter = 1
+ for i in eachindex(model.reactions)
+ if i in objective_indices
+ opt_weights[i] = weights[wcounter]
+ wcounter += 1
+ end
+ end
+
+ v = opt_model[:x]
+ @objective(opt_model, sense, sum(opt_weights[i] * v[i] for i in objective_indices))
+ end
+end
+
+"""
+ modify_solver(optimizer)
+
+Modify the solver used to solve the model.
+Typically the solver is specified as a required argument.
+This function is useful if the problem has multiple subparts that require different solvers.
+See
+"""
+function modify_solver(optimizer)
+ (model, opt_model) -> begin
+ JuMP.set_optimizer(opt_model, optimizer)
+ end
+end
diff --git a/src/base/solver.jl b/src/base/solver.jl
index fd7b6ab7dc8994e0f4fde6eeab92f9e9fecff12a..0a159d622cfafecbf961c5feeab02fea49a56778 100644
--- a/src/base/solver.jl
+++ b/src/base/solver.jl
@@ -29,7 +29,6 @@ function optimize_model(
sense = MOI.MIN_SENSE,
) where {LM<:MetabolicModel}
optimization_model = make_optimization_model(model, optimizer; sense = sense)
- x = optimization_model[:x]
JuMP.optimize!(optimization_model)
return optimization_model
end
diff --git a/src/base/utilities.jl b/src/base/utilities.jl
index e2d6d7233bd8234886c830f7f03c94f03f50bb34..59f202da7e1c259e25bac84fdb206e92d8e3b528 100644
--- a/src/base/utilities.jl
+++ b/src/base/utilities.jl
@@ -18,31 +18,3 @@ Base.isequal(model1::CoupledLinearModel, model2::CoupledLinearModel) =
Base.copy(model::CoupledLinearModel) =
CoupledLinearModel(model.lm, model.C, model.cl, model.cu)
-
-
-"""
- set_bound(index, optimization_model; ub=1000, lb=-1000)
-Helper function to set the bounds of variables.
-The JuMP `set_normalized_rhs` function is a little confusing, so this function simplifies setting
-constraints. Just supply the constraint `index` and the model and they will be set to `ub` and `lb`.
-"""
-function set_bound(vind, opt_model; ub = 1000, lb = -1000)
- if lb <= 0
- set_normalized_rhs(opt_model[:lbs][vind], abs(lb))
- else
- set_normalized_rhs(opt_model[:lbs][vind], -abs(lb))
- end
- set_normalized_rhs(opt_model[:ubs][vind], ub)
-end
-
-"""
- modify_constraint(reaction::Reaction, lb, ub)
-
-Modify constraints of model reaction.
-"""
-function modify_constraint(reaction::Reaction, lb, ub)
- (model, opt_model) -> begin
- ind = model.reactions[reaction]
- set_bound(ind, opt_model, lb=lb, ub=ub)
- end
-end
diff --git a/src/io/io.jl b/src/io/io.jl
index 843d42ddee4001bc454fd90f4a01ac9630b195ef..ba12321a4848521db08670fdcbb61f6b14b65a37 100644
--- a/src/io/io.jl
+++ b/src/io/io.jl
@@ -1,7 +1,7 @@
"""
read_model(file_location::String))
-Reads a model at `file_location` and returns a constraint based `model::CobraModel`.
+Reads a model at `file_location` and returns a constraint based `model::StandardModel`.
Currently supported formats include SBML (.xml), Matlab (.mat) and JSON (.json) models.
The model format is inferred from the `file_location` extension.
@@ -25,31 +25,31 @@ function read_model(file_location::String)
model = read_json_model(file_location)
catch err
@error "JSON model reading error.\n$err"
- model = CobraModel()
+ model = StandardModel()
end
elseif endswith(file_location, ".xml")
try
model = reconstruct_model_sbml(file_location)
catch err
@error "SBML model reading error.\n$err"
- model = CobraModel()
+ model = StandardModel()
end
elseif endswith(file_location, ".mat")
try
model = read_matlab_model(file_location)
catch err
@error "Matlab model reading error.\n$err"
- model = CobraModel()
+ model = StandardModel()
end
else
@error "Model format not supported. The format is inferred from the file extension. Supported formats: *.mat, *.xml, *.json."
- model = CobraModel()
+ model = StandardModel()
end
return model
end
"""
- save_model(model::CobraModel, file_location::String)
+ save_model(model::StandardModel, file_location::String)
Save model at `file_location`. Infers format from `file_location` extension.
Supported formats include SBML (.xml), Matlab COBRA models (.mat) and JSON COBRA models (.json).
@@ -57,7 +57,7 @@ Supported formats include SBML (.xml), Matlab COBRA models (.mat) and JSON COBRA
Note, only the fields contained in model are saved. Make sure that information isn't
lost between reading a model and writing a model (e.g. check gene reaction rules, notes and annotations).
"""
-function save_model(model::CobraModel, file_location::String)
+function save_model(model::StandardModel, file_location::String)
if endswith(file_location, ".json")
write_json_model(model, file_location)
elseif endswith(file_location, ".xml")
@@ -126,13 +126,13 @@ function reconstruct_model_sbml(file_location::String)
# m.species
# m.compartments
# return Model()
- return CobraModel()
+ return StandardModel()
end
"""
- save_sbml_model(model::CobraModel, file_location::String)
+ save_sbml_model(model::StandardModel, file_location::String)
"""
-function save_sbml_model(model::CobraModel, file_location::String)
+function save_sbml_model(model::StandardModel, file_location::String)
# To do...
end
diff --git a/src/io/json_reader.jl b/src/io/json_reader.jl
index c767b97917c74235e99209d54d8f3160697d99aa..222604d363b8a9bb0f7ce0784d33ea4adb17c6e2 100644
--- a/src/io/json_reader.jl
+++ b/src/io/json_reader.jl
@@ -138,5 +138,5 @@ function read_json_model(file_location::String)
)
end
- return CobraModel(modelid, rxns, mets, genes)
+ return StandardModel(modelid, rxns, mets, genes)
end
diff --git a/src/io/json_writer.jl b/src/io/json_writer.jl
index a0686de759348e54a0b2b1180d22c2adf6d719e5..3aac7c956fff13566df4cad3e8ac5e82827603a6 100644
--- a/src/io/json_writer.jl
+++ b/src/io/json_writer.jl
@@ -1,7 +1,7 @@
"""
- save_json_model(model::CobraModel, file_location::String)
+ save_json_model(model::StandardModel, file_location::String)
"""
-function write_json_model(model::CobraModel, file_location::String)
+function write_json_model(model::StandardModel, file_location::String)
modeldict = Dict{String,Any}()
modeldict["id"] = model.id
diff --git a/src/io/m_reader.jl b/src/io/m_reader.jl
index d1427c50b26a5e71f28dc305a9bc35bff1d7a8cd..8c5379244c4fbf737529cec5a1b63c805ca3eb77 100644
--- a/src/io/m_reader.jl
+++ b/src/io/m_reader.jl
@@ -136,5 +136,5 @@ function read_matlab_model(file_location::String)
)
end
- return CobraModel(model_id, rxns, mets, genes)
+ return StandardModel(model_id, rxns, mets, genes)
end
diff --git a/src/io/m_writer.jl b/src/io/m_writer.jl
index b646e1211d0aae1d07e9d3fa523b97ea68baf985..3938bbc797ab24f4a1cdb976b51e7e3af3cc447d 100644
--- a/src/io/m_writer.jl
+++ b/src/io/m_writer.jl
@@ -1,9 +1,9 @@
"""
- write_matlab_model(model::CobraModel, file_location::String)
+ write_matlab_model(model::StandardModel, file_location::String)
Some information is lost here, e.g. notes and some annotations.
"""
-function write_matlab_model(model::CobraModel, file_location::String)
+function write_matlab_model(model::StandardModel, file_location::String)
S = stoichiometry(model)
b = balance(model)
lbs, ubs = bounds(model)
diff --git a/src/reconstruction/model_manipulations.jl b/src/reconstruction/model_manipulations.jl
index d6b6c3ed9abe128e09840de5995838e904441329..4aa977ef4d57f83deb69c3476b7eefe0171d3678 100644
--- a/src/reconstruction/model_manipulations.jl
+++ b/src/reconstruction/model_manipulations.jl
@@ -1,15 +1,15 @@
"""
- add!(model::CobraModel, rxns::Union{Array{Reaction, 1}, Reaction})
+ add!(model::StandardModel, rxns::Union{Array{Reaction, 1}, Reaction})
Add `rxn(s)` to `model` if they are not already present based on reaction `id`.
"""
-function add!(model::CobraModel, rxns::Array{Reaction,1})
+function add!(model::StandardModel, rxns::Array{Reaction,1})
for rxn in rxns
add!(model, rxn)
end
end
-function add!(model::CobraModel, rxn::Reaction)
+function add!(model::StandardModel, rxn::Reaction)
if model[rxn] == -1
push!(model.reactions, rxn)
else
@@ -19,17 +19,17 @@ function add!(model::CobraModel, rxn::Reaction)
end
"""
- add!(model::CobraModel, mets::Union{Array{Metabolite, 1}, Metabolite})
+ add!(model::StandardModel, mets::Union{Array{Metabolite, 1}, Metabolite})
Add `met(s)` to `model` if they are not already present, based on metabolite `id`.
"""
-function add!(model::CobraModel, mets::Array{Metabolite,1})
+function add!(model::StandardModel, mets::Array{Metabolite,1})
for met in mets
add!(model, met)
end
end
-function add!(model::CobraModel, met::Metabolite)
+function add!(model::StandardModel, met::Metabolite)
if model[met] == -1
push!(model.metabolites, met)
else
@@ -39,17 +39,17 @@ function add!(model::CobraModel, met::Metabolite)
end
"""
- add!(model::CobraModel, genes::Union{Array{Gene, 1}, Gene})
+ add!(model::StandardModel, genes::Union{Array{Gene, 1}, Gene})
Add `gene(s)` to `model` if they are not already present based on gene `id`.
"""
-function add!(model::CobraModel, genes::Array{Gene,1})
+function add!(model::StandardModel, genes::Array{Gene,1})
for gene in genes
add!(model, gene)
end
end
-function add!(model::CobraModel, gene::Gene)
+function add!(model::StandardModel, gene::Gene)
if model[gene] == -1
push!(model.genes, gene)
else
@@ -59,11 +59,11 @@ function add!(model::CobraModel, gene::Gene)
end
"""
- rm!(model::CobraModel, rxns::Union{Array{Reaction, 1}, Reaction})
+ rm!(model::StandardModel, rxns::Union{Array{Reaction, 1}, Reaction})
Remove all `rxn(s)` from `model` if the `id`s match those in `rxns`.
"""
-function rm!(model::CobraModel, rxns::Union{Array{Reaction,1},Reaction})
+function rm!(model::StandardModel, rxns::Union{Array{Reaction,1},Reaction})
new_rxn_list = Reaction[]
for r in model.reactions
if typeof(rxns) == Reaction
@@ -81,11 +81,11 @@ function rm!(model::CobraModel, rxns::Union{Array{Reaction,1},Reaction})
end
"""
- rm!(model::CobraModel, mets::Union{Array{Metabolite, 1}, Metabolite})
+ rm!(model::StandardModel, mets::Union{Array{Metabolite, 1}, Metabolite})
Remove `met(s)` from `model` based on metabolite `id`.
"""
-function rm!(model::CobraModel, mets::Union{Array{Metabolite,1},Metabolite})
+function rm!(model::StandardModel, mets::Union{Array{Metabolite,1},Metabolite})
new_met_list = Metabolite[]
for m in model.metabolites
if typeof(mets) == Metabolite
@@ -103,11 +103,11 @@ function rm!(model::CobraModel, mets::Union{Array{Metabolite,1},Metabolite})
end
"""
- rm!(model::CobraModel, genes::Union{Array{Gene, 1}, Gene})
+ rm!(model::StandardModel, genes::Union{Array{Gene, 1}, Gene})
Remove `gene(s)` from `model` based on gene `id`.
"""
-function rm!(model::CobraModel, genes::Union{Array{Gene,1},Gene})
+function rm!(model::StandardModel, genes::Union{Array{Gene,1},Gene})
new_gene_list = Gene[]
for g in model.genes
if typeof(genes) == Gene
@@ -125,7 +125,7 @@ function rm!(model::CobraModel, genes::Union{Array{Gene,1},Gene})
end
"""
- fix_model!(model::CobraModel)
+ fix_model!(model::StandardModel)
Inspect metabolites and genes of `model` relative to the reactions of `model`.
Remove genes or metabolites that are not used in the reactions.
@@ -135,7 +135,7 @@ a different `id`. Other functions are provided to help identify these cases.
See also: [`check_duplicate_annotations`](@ref), [`check_duplicate_reaction`](@ref), [`check_same_formula`](@ref).
"""
-function fix_model!(model::CobraModel)
+function fix_model!(model::StandardModel)
rxn_mets = Metabolite[] # list of metabolites used in reactions
for rxn in model.reactions
for met in keys(rxn.metabolites)
diff --git a/src/sampling/sampling_tools.jl b/src/sampling/sampling_tools.jl
index 90551e75f37a9118c806d39a919d72649bee5e61..991965565e2951ad72805b155da469250ad64e7b 100644
--- a/src/sampling/sampling_tools.jl
+++ b/src/sampling/sampling_tools.jl
@@ -3,17 +3,21 @@
Generate warmup points for all the reactions in the model that
are not fixed. Assumes you feed in a JuMP model that is already
-constrained i.e. the constrains are already applied into cbm.
-Note, extra constraints applied to ubs and lbs will have no effect.
+constrained i.e. the constraints are already applied to `cbm`.
-numstop (maximum value is 2*number of warmup points). Reduce this to prematurely
-stop finding warmup points.
+If you do not want to generate all possible warmup points
+(2*n where n is the number of non-fixed variables), then reduce `numstop`
+to the desired number of warmup points. This number of warmup points will
+be randomly drawn from the non-fixed fluxes.
+
+By default each flux will be minimized and maximized, however, random objectives
+can be used by setting `random_objective` true.
"""
function get_warmup_points(cbm; random_objective = false, numstop = 1e10)
v = cbm[:x]
ubs = cbm[:ubs]
lbs = cbm[:lbs]
- # determine which rxns should be max/min-ized
+ # determine which rxns should be max/min-ized (non fixed fluxes)
fixed_rxns = Int64[]
for i in eachindex(v)
ub_val = normalized_rhs(ubs[i])
@@ -25,11 +29,12 @@ function get_warmup_points(cbm; random_objective = false, numstop = 1e10)
end
end
- # determine number of warmup points
+ # determine number of warmup points (if numstops is less than the number of non-fixed reactions.)
var_rxn_inds = shuffle!(filter(x -> !(x in fixed_rxns), 1:length(v))) # shuffle the optimization points
NN = numstop > length(var_rxn_inds) ? length(var_rxn_inds) : numstop
wpoints = zeros(length(v), 2 * NN)
+ # generate warmup points.
for (i, ii) in zip(1:length(var_rxn_inds), 1:2:(2*length(var_rxn_inds)))
i > NN && break
@@ -62,7 +67,8 @@ end
"""
get_bound_vectors(opt_model)
-Return Float64 vectors of the lower and upper bounds of the JuMP constraint refs.
+Return Float64 vectors of the lower and upper bounds of `opt_model` constraints,
+which is a JuMP model.
"""
function get_bound_vectors(opt_model)
lbconref = opt_model[:lbs]
@@ -77,17 +83,19 @@ function get_bound_vectors(opt_model)
end
end
ubs = [normalized_rhs(ubconref[i]) for i in eachindex(ubconref)]
+
return lbs, ubs
end
"""
- hit_and_run(N::Int64, model::CobraModel, optimizer; constraints=Dict{String, Tuple{Float64,Float64}}(), keepevery=100, samplesize=1000, solver_attributes=Dict{Any, Any}(), random_objective=false)
+ hit_and_run(N::Int64, cbm; keepevery=100, samplesize=1000, random_objective=false)
-Perform basic hit and run sampling for `N` iterations using `model` with `optimizer` from `JuMP`.
-Additional constraints supplied by `constraints` as a dictionary of reaction `id`s mapped to a tuple of `(lb, ub)` of fluxes.
+Perform basic hit and run sampling for `N` iterations on `cbm`.
+Note that `cbm` is a JuMP model that contains whatever solver will be used to find the warmup points.
+Also note, `cbm` should already be fully constrained by however you desire.
Every `keepevery` iteration is logged as a sample, where the sample size matrix has `samplesize` columns.
-Solver specific settings can be set using `solver_attributes`.
-Warm up points are generated in a flux variability sense, unless `random_objective` is true, in which case a randomly weighted objective is used 2*number of reactions to define the warmup points.
+Warm up points are generated in a flux variability sense, unless `random_objective` is true,
+in which case a randomly weighted objective is used 2*number of reactions to define the warmup points.
Note that N needs to be >> samplesize.
Sample size is the size of the samples kept in memory.
@@ -97,31 +105,11 @@ See also: [`achr`](@ref)
"""
function hit_and_run(
N::Int64,
- model::CobraModel,
- optimizer;
- constraints = Dict{String,Tuple{Float64,Float64}}(),
+ cbm;
keepevery = 100,
samplesize = 1000,
- solver_attributes = Dict{Any,Any}(),
random_objective = false,
- sense = MOI.MAX_SENSE,
)
- # get core optimization problem
- 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
- 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
-
lbs, ubs = get_bound_vectors(cbm) # get actual ub and lb constraints, can't use model function because the user may have changed them in the function arguments
wpoints = get_warmup_points(cbm; random_objective = random_objective)
@@ -185,25 +173,19 @@ function hit_and_run(
end
- violation_inds = test_samples(samples, model, ubs, lbs)
- if !isempty(violation_inds)
- @warn "Samples: $violation_inds do not satisfy the problem constraints."
- end
-
return samples
end
"""
- test_samples(samples::Array{Float64, 2}, model::CobraModel, ubs, lbs)
+ test_samples(samples::Array{Float64, 2}, model::StandardModel, ubs, lbs)
Test if samples pass tests: mass balances and constraints are satisfied..
"""
-function test_samples(samples::Array{Float64,2}, model::CobraModel, ubs, lbs)
- S = stoichiometry(model) # assume S has not been modified from model AND b is zero
+function test_samples(samples::Array{Float64,2}, mass_balance, balance, lbs, ubs)
violations = Int64[]
tol = 1e-6
for i = 1:size(samples, 2)
- if sum(abs, S * samples[:, i]) < tol
+ if sum(abs, mass_balance * samples[:, i] - balance) < tol
equality = true
else
equality = false
@@ -219,106 +201,106 @@ function test_samples(samples::Array{Float64,2}, model::CobraModel, ubs, lbs)
return violations
end
-"""
- achr(N::Int64, model::CobraModel, optimizer;constraints=Dict{String, Tuple{Float64,Float64}}(), keepevery=100, samplesize=1000, solver_attributes=Dict{Any, Any}(), random_objective=false)
-
-Perform artificially centered hit and run.
-Uses the same arguments as the `hit_and_run` sampler.
-Needs work, for long iterations it becomes unstable (violates bounds).
-
-See also: [`hit_and_run`](@ref)
-"""
-function achr(
- N::Int64,
- model::CobraModel,
- optimizer;
- constraints = Dict{String,Tuple{Float64,Float64}}(),
- keepevery = 100,
- samplesize = 1000,
- solver_attributes = Dict{Any,Any}(),
- random_objective = false,
-)
- cbm, v, _, ubcons, lbcons = build_cbm(model)
-
- set_optimizer(cbm, optimizer) # choose optimizer
- 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, ubcons, lbcons; ub = con[1], lb = con[2])
- end
-
- ubs, lbs = get_bound_vectors(ubcons, lbcons) # get actual ub and lb constraints
-
- wpoints = get_warmup_points(cbm, v, ubcons, lbcons, random_objective = random_objective)
-
- nwpts = size(wpoints, 2) # number of warmup points generated
- samples = zeros(size(wpoints, 1), samplesize) # sample storage
- current_point = zeros(size(wpoints, 1))
- current_point .= wpoints[:, rand(1:nwpts)] # pick random initial point
-
- shat = mean(wpoints, dims = 2)[:] # mean point
-
- δdirtol = 1e-6 # too small directions get ignored ≈ 0 (solver precision issue)
- sample_num = 0
- samplelength = 0
- updatesamplesizelength = true
- for n = 1:N
- if updatesamplesizelength # switch to samples
- direction_point = (@view wpoints[:, rand(1:nwpts)]) - (@view current_point[:]) # use warmup points to find direction in warmup phase
- else
- direction_point = (@view samples[:, rand(1:(samplelength))]) - (@view shat[:]) # after warmup phase, only find directions in sampled space
- end
-
- λmax = 1e10
- λmin = -1e10
- for i in eachindex(lbs)
- δlower = lbs[i] - current_point[i]
- δupper = ubs[i] - current_point[i]
- if direction_point[i] < -δdirtol
- lower = δupper / direction_point[i]
- upper = δlower / direction_point[i]
- elseif direction_point[i] > δdirtol
- lower = δlower / direction_point[i]
- upper = δupper / direction_point[i]
- else
- lower = -1e10
- upper = 1e10
- end
- lower > λmin && (λmin = lower) # max min step size that satisfies all bounds
- upper < λmax && (λmax = upper) # min max step size that satisfies all bounds
- end
-
- if λmax <= λmin || λmin == -1e10 || λmax == 1e10 # this sometimes can happen
- @warn "Infeasible direction at iteration $(n)..."
- continue
- end
-
- λ = rand() * (λmax - λmin) + λmin # random step size
- current_point .= current_point .+ λ .* direction_point # will be feasible
- shat .= (n .* shat .+ current_point) ./ (n + 1)
-
- if n % keepevery == 0
- sample_num += 1
- samples[:, sample_num] .= current_point
- if sample_num >= samplesize
- updatesamplesizelength = false
- sample_num = 0 # reset, start replacing the older samples
- end
- updatesamplesizelength && (samplelength += 1) # lags sample_num because the latter is a flag as well
- end
-
- end
-
- violation_inds = test_samples(samples, model, ubs, lbs)
- if !isempty(violation_inds)
- @warn "Samples: $violation_inds do not satisfy the problem constraints."
- end
-
- return samples
-end
+# """
+# achr(N::Int64, model::StandardModel, optimizer;constraints=Dict{String, Tuple{Float64,Float64}}(), keepevery=100, samplesize=1000, solver_attributes=Dict{Any, Any}(), random_objective=false)
+
+# Perform artificially centered hit and run.
+# Uses the same arguments as the `hit_and_run` sampler.
+# Needs work, for long iterations it becomes unstable (violates bounds).
+
+# See also: [`hit_and_run`](@ref)
+# """
+# function achr(
+# N::Int64,
+# model::StandardModel,
+# optimizer;
+# constraints = Dict{String,Tuple{Float64,Float64}}(),
+# keepevery = 100,
+# samplesize = 1000,
+# solver_attributes = Dict{Any,Any}(),
+# random_objective = false,
+# )
+# cbm, v, _, ubcons, lbcons = build_cbm(model)
+
+# set_optimizer(cbm, optimizer) # choose optimizer
+# 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, ubcons, lbcons; ub = con[1], lb = con[2])
+# end
+
+# ubs, lbs = get_bound_vectors(ubcons, lbcons) # get actual ub and lb constraints
+
+# wpoints = get_warmup_points(cbm, v, ubcons, lbcons, random_objective = random_objective)
+
+# nwpts = size(wpoints, 2) # number of warmup points generated
+# samples = zeros(size(wpoints, 1), samplesize) # sample storage
+# current_point = zeros(size(wpoints, 1))
+# current_point .= wpoints[:, rand(1:nwpts)] # pick random initial point
+
+# shat = mean(wpoints, dims = 2)[:] # mean point
+
+# δdirtol = 1e-6 # too small directions get ignored ≈ 0 (solver precision issue)
+# sample_num = 0
+# samplelength = 0
+# updatesamplesizelength = true
+# for n = 1:N
+# if updatesamplesizelength # switch to samples
+# direction_point = (@view wpoints[:, rand(1:nwpts)]) - (@view current_point[:]) # use warmup points to find direction in warmup phase
+# else
+# direction_point = (@view samples[:, rand(1:(samplelength))]) - (@view shat[:]) # after warmup phase, only find directions in sampled space
+# end
+
+# λmax = 1e10
+# λmin = -1e10
+# for i in eachindex(lbs)
+# δlower = lbs[i] - current_point[i]
+# δupper = ubs[i] - current_point[i]
+# if direction_point[i] < -δdirtol
+# lower = δupper / direction_point[i]
+# upper = δlower / direction_point[i]
+# elseif direction_point[i] > δdirtol
+# lower = δlower / direction_point[i]
+# upper = δupper / direction_point[i]
+# else
+# lower = -1e10
+# upper = 1e10
+# end
+# lower > λmin && (λmin = lower) # max min step size that satisfies all bounds
+# upper < λmax && (λmax = upper) # min max step size that satisfies all bounds
+# end
+
+# if λmax <= λmin || λmin == -1e10 || λmax == 1e10 # this sometimes can happen
+# @warn "Infeasible direction at iteration $(n)..."
+# continue
+# end
+
+# λ = rand() * (λmax - λmin) + λmin # random step size
+# current_point .= current_point .+ λ .* direction_point # will be feasible
+# shat .= (n .* shat .+ current_point) ./ (n + 1)
+
+# if n % keepevery == 0
+# sample_num += 1
+# samples[:, sample_num] .= current_point
+# if sample_num >= samplesize
+# updatesamplesizelength = false
+# sample_num = 0 # reset, start replacing the older samples
+# end
+# updatesamplesizelength && (samplelength += 1) # lags sample_num because the latter is a flag as well
+# end
+
+# end
+
+# violation_inds = test_samples(samples, model, ubs, lbs)
+# if !isempty(violation_inds)
+# @warn "Samples: $violation_inds do not satisfy the problem constraints."
+# end
+
+# return samples
+# end
diff --git a/src/types/agene.jl b/src/types/gene.jl
similarity index 100%
rename from src/types/agene.jl
rename to src/types/gene.jl
diff --git a/src/types/ametabolite.jl b/src/types/metabolite.jl
similarity index 100%
rename from src/types/ametabolite.jl
rename to src/types/metabolite.jl
diff --git a/src/types/areaction.jl b/src/types/reaction.jl
similarity index 100%
rename from src/types/areaction.jl
rename to src/types/reaction.jl
diff --git a/src/types/metabolicModel.jl b/src/types/standardModel.jl
similarity index 68%
rename from src/types/metabolicModel.jl
rename to src/types/standardModel.jl
index f2c354d418f0e5354ed9d2662fde076127e26488..36a70c14b27f7b461e970e2a40df4d71836e4da8 100644
--- a/src/types/metabolicModel.jl
+++ b/src/types/standardModel.jl
@@ -1,5 +1,5 @@
"""
-CobraModel struct of a constraint based metabolic model.
+StandardModel struct of a constraint based metabolic model.
# Fields
````
@@ -9,7 +9,7 @@ metabolites :: Array{Metabolite, 1}
genes :: Array{Gene, 1}
````
"""
-mutable struct CobraModel <: MetabolicModel
+mutable struct StandardModel <: MetabolicModel
id::String
reactions::Array{Reaction,1}
metabolites::Array{Metabolite,1}
@@ -17,9 +17,9 @@ mutable struct CobraModel <: MetabolicModel
end
"""
-Pretty printing of model::CobraModel.
+Pretty printing of model::StandardModel.
"""
-function Base.show(io::IO, ::MIME"text/plain", m::CobraModel)
+function Base.show(io::IO, ::MIME"text/plain", m::StandardModel)
println(
io,
"Constraint based model: ",
@@ -37,61 +37,61 @@ function Base.show(io::IO, ::MIME"text/plain", m::CobraModel)
end
"""
-CobraModel()
+StandardModel()
Empty model constructor.
"""
-function CobraModel()
- CobraModel("blank", Array{Reaction,1}(), Array{Metabolite,1}(), Array{Gene,1}())
+function StandardModel()
+ StandardModel("blank", Array{Reaction,1}(), Array{Metabolite,1}(), Array{Gene,1}())
end
"""
- getindex(model::CobraModel, rxn::Reaction)
+ getindex(model::StandardModel, rxn::Reaction)
Get the index of `rxn` in `model`, based on reaction `id`.
Return -1 if not found.
Typical usage: ind = model[rxn]
"""
-function Base.getindex(model::CobraModel, rxn::Reaction)
+function Base.getindex(model::StandardModel, rxn::Reaction)
return model.reactions[rxn]
end
"""
- getindex(model::CobraModel, met::Metabolite)
+ getindex(model::StandardModel, met::Metabolite)
Get the index of `met` in `model`, based on metabolite `id`.
Return -1 if not found.
Typical usage: ind = model[met]
"""
-function Base.getindex(model::CobraModel, met::Metabolite)
+function Base.getindex(model::StandardModel, met::Metabolite)
return model.metabolites[met]
end
"""
- getindex(model::CobraModel, gene::Gene)
+ getindex(model::StandardModel, gene::Gene)
Get the index of `gene` in `model`, based on gene `id`.
Return -1 if not found.
Typical usage: ind = model[gene]
"""
-function Base.getindex(model::CobraModel, gene::Gene)
+function Base.getindex(model::StandardModel, gene::Gene)
return model.genes[gene]
end
# generic interface functions
-function reactions(model::CobraModel)::Vector{String}
+function reactions(model::StandardModel)::Vector{String}
[r.id for r in model.reactions]
end
-function metabolites(model::CobraModel)::Vector{String}
+function metabolites(model::StandardModel)::Vector{String}
[m.id for m in model.metabolites]
end
-function stoichiometry(model::CobraModel)::SparseMat
+function stoichiometry(model::StandardModel)::SparseMat
S = SparseArrays.spzeros(length(model.metabolites), length(model.reactions))
metids = metabolites(model)
for (i, rxn) in enumerate(model.reactions) # column
@@ -106,17 +106,17 @@ function stoichiometry(model::CobraModel)::SparseMat
return S
end
-function bounds(model::CobraModel)::Tuple{SparseVec,SparseVec}
+function bounds(model::StandardModel)::Tuple{SparseVec,SparseVec}
ubs = [rxn.ub for rxn in model.reactions]
lbs = [rxn.lb for rxn in model.reactions]
return lbs, ubs
end
-function balance(model::CobraModel)::SparseVec
+function balance(model::StandardModel)::SparseVec
SparseArrays.spzeros(length(model.metabolites))
end
-function objective(model::CobraModel)::SparseVec
+function objective(model::StandardModel)::SparseVec
obj_arr = SparseArrays.spzeros(length(model.reactions))
j = -1
for (i, r) in enumerate(model.reactions)
diff --git a/test/analysis/fba.jl b/test/analysis/fba.jl
index 4e89af7b5cf72bf4813216002029669f6927571c..86fa3bca8fd79746e9e223f874b23960959ae18e 100644
--- a/test/analysis/fba.jl
+++ b/test/analysis/fba.jl
@@ -1,23 +1,20 @@
-@testset "Flux balance analysis" begin
+@testset "Flux balance analysis with LinearModel" begin
cp = test_simpleLP()
lp = flux_balance_analysis(cp, GLPK.Optimizer)
- x = lp[:x]
@test termination_status(lp) === MOI.OPTIMAL
- sol = COBREXA.JuMP.value.(x)
+ sol = JuMP.value.(lp[:x])
@test sol ≈ [1.0, 2.0]
lp = flux_balance_analysis(cp, Clp.Optimizer)
- x = lp[:x]
@test termination_status(lp) === MOI.OPTIMAL
- sol = COBREXA.JuMP.value.(x)
+ sol = JuMP.value.(lp[:x])
@test sol ≈ [1.0, 2.0]
# test the maximization of the objective
cp = test_simpleLP2()
lp = flux_balance_analysis(cp, GLPK.Optimizer)
- x = lp[:x]
@test termination_status(lp) === MOI.OPTIMAL
- sol = COBREXA.JuMP.value.(x)
+ sol = JuMP.value.(lp[:x])
@test sol ≈ [-1.0, 2.0]
# test with a more biologically meaningfull model
@@ -31,9 +28,8 @@
expected_optimum = 0.9219480950504393
lp = flux_balance_analysis(cp, GLPK.Optimizer)
- x = lp[:x]
@test termination_status(lp) === MOI.OPTIMAL
- sol = COBREXA.JuMP.value.(x)
+ sol = JuMP.value.(lp[:x])
@test objective_value(lp) ≈ expected_optimum
@test cp.c' * sol ≈ expected_optimum
@@ -45,7 +41,7 @@
@test all([fluxes_dict[rxns[i]] == sol[i] for i in eachindex(rxns)])
end
-@testset "Flux balance analysis with CobraModel" begin
+@testset "Flux balance analysis with StandardModel" begin
model = read_model(
download_data_file(
"http://bigg.ucsd.edu/static/models/e_coli_core.json",
@@ -53,49 +49,15 @@ end
"7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
),
)
- @test length(model.reactions) == 95 # read in correctly
- # FBA
biomass = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")
- cons = Dict("EX_glc__D_e" => (-12.0, -12.0))
- optimizer = COBREXA.Tulip.Optimizer # quiet by default
- sol = fba(model, optimizer; objective_func = biomass, constraints = cons)
- pfl = findfirst(model.reactions, "PFL")
- solmulti = fba(model, optimizer; objective_func = [biomass, pfl], weights = [0.8, 0.2]) # classic flux balance analysis
-
-
- flux_vec = [sol[rxn.id] for rxn in model.reactions]
- sol_mapped = map_fluxes(flux_vec, model)
- @test isapprox(sol_mapped["BIOMASS_Ecoli_core_w_GAM"], 1.0572509997013568, atol = 1e-6)
+ glucose = findfirst(model.reactions, "EX_glc__D_e")
+ sol = flux_balance_analysis_dict(model, Tulip.Optimizer; modifications=[modify_objective(biomass), modify_constraint(glucose, -12, -12), modify_sense(MOI.MAX_SENSE), modify_solver_attribute("IPM_IterationsLimit", 110)])
@test isapprox(sol["BIOMASS_Ecoli_core_w_GAM"], 1.0572509997013568, atol = 1e-6)
- @test !isempty(solmulti)
-
- sol = fba(model, optimizer; objective_func = biomass)
- # atom tracker
- atom_fluxes = atom_exchange(sol, model)
- @test isapprox(atom_fluxes["C"], -37.1902, atol = 1e-3)
-
- # exchange trackers
- consuming, producing = exchange_reactions(sol; verbose = false)
- @test isapprox(consuming["EX_nh4_e"], -4.76532, atol = 1e-3)
-
- # metabolite trackers
- consuming, producing = metabolite_fluxes(sol, model)
- @test isapprox(consuming["atp_c"]["PFK"], -7.47738, atol = 1e-3)
- @test isapprox(producing["atp_c"]["PYK"], 1.75818, atol = 1e-3)
-
- # set bounds
- cbm = make_optimization_model(model, optimizer)
- ubs = cbm[:ubs]
- lbs = cbm[:lbs]
- glucose_index = model[findfirst(model.reactions, "EX_glc__D_e")]
- o2_index = model[findfirst(model.reactions, "EX_o2_e")]
- atpm_index = model[findfirst(model.reactions, "ATPM")]
- set_bound(glucose_index, cbm; ub = -1.0, lb = -1.0)
- @test normalized_rhs(ubs[glucose_index]) == -1.0
- @test normalized_rhs(lbs[glucose_index]) == 1.0
- set_bound(o2_index, cbm; ub = 1.0, lb = 1.0)
- @test normalized_rhs(ubs[o2_index]) == 1.0
- @test normalized_rhs(lbs[o2_index]) == -1.0
+ pfl = findfirst(model.reactions, "PFL")
+ pfl_frac = 0.8
+ biomass_frac = 0.2
+ sol_multi = flux_balance_analysis_dict(model, Tulip.Optimizer; modifications= modify_objective([biomass, pfl]; weights = [biomass_frac, pfl_frac]))
+ @test isapprox(biomass_frac*sol_multi["BIOMASS_Ecoli_core_w_GAM"] + pfl_frac*sol_multi["PFL"], 31.999999998962604, atol = 1e-6)
end
diff --git a/test/analysis/fva.jl b/test/analysis/fva.jl
index 8c555dc8c2482d6ba3377a5c7396f81b23c68b80..17fa868ad3c249042b10cb8aebb866415afb3719 100644
--- a/test/analysis/fva.jl
+++ b/test/analysis/fva.jl
@@ -60,7 +60,7 @@ end
rmprocs(pids)
end
-@testset "Flux variability analysis with CobraModel" begin
+@testset "Flux variability analysis with StandardModel" begin
model = read_model(
download_data_file(
"http://bigg.ucsd.edu/static/models/e_coli_core.json",
@@ -68,28 +68,10 @@ end
"7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
),
)
- @test length(model.reactions) == 95 # read in correctly
biomass = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")
- pfl = findfirst(model.reactions, "PFL")
-
- # FVA
- optimizer = COBREXA.Tulip.Optimizer
- atts = Dict("IPM_IterationsLimit" => 500)
- cons = Dict("EX_glc__D_e" => (-10.0, -10.0))
- fva_max, fva_min =
- fva(model, optimizer; objective_func = biomass, solver_attributes = atts)
- fva_max2, fva_min2 = fva(
- model,
- optimizer;
- objective_func = [biomass, pfl],
- 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
+ glucose = findfirst(model.reactions, "EX_glc__D_e")
+ fva_max, fva_min = flux_variability_analysis(model, Tulip.Optimizer; modifications=[modify_solver_attribute("IPM_IterationsLimit", 500), modify_constraint(glucose, -10, -10)])
+ @test isapprox(fva_max["PDH"]["PDH"], 9.338920475333763, atol = 1e-6)
+ @test isapprox(fva_min["PDH"]["PDH"], 9.270274952732315, atol = 1e-6)
end
diff --git a/test/analysis/pfba.jl b/test/analysis/pfba.jl
index fffd674fe770ba92855c70b2b18860e5336ac53b..106e5aa6512f82a8bf5621c8d1e21e7eb1454ffe 100644
--- a/test/analysis/pfba.jl
+++ b/test/analysis/pfba.jl
@@ -1,5 +1,4 @@
-@testset "Parsimonious flux balance analysis with CobraModel" begin
-
+@testset "Parsimonious flux balance analysis with StandardModel" begin
model = read_model(
download_data_file(
"http://bigg.ucsd.edu/static/models/e_coli_core.json",
@@ -9,29 +8,11 @@
)
biomass = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")
- cons = Dict("EX_glc__D_e" => (-12.0, -12.0))
-
- # pFBA
- atts = Dict(
- "eps_abs" => 5e-4,
- "eps_rel" => 5e-4,
- "max_iter" => 100_000,
- "verbose" => false,
- )
- solworks = pfba(
- model,
- OSQP.Optimizer;
- objective_func = biomass,
- solver_attributes = atts,
- constraints = cons,
- ) # just see if it works - OSQP is a terrible LP solver
- sol = pfba(
- model,
- [COBREXA.Tulip.Optimizer, COBREXA.OSQP.Optimizer];
- objective_func = biomass,
- solver_attributes = Dict("opt1" => Dict{Any,Any}(), "opt2" => atts),
- ) # try two optimizers
-
- @test !isempty(solworks)
- @test isapprox(sol["PGM"], -14.737442319041387, atol = 1e-3)
+ glucose = findfirst(model.reactions, "EX_glc__D_e")
+
+ d = parsimonious_flux_balance_analysis_dict(model, Tulip.Optimizer; modifications=modify_constraint(glucose, -12, -12), qp_solver=modify_solver(OSQP.Optimizer), qp_solver_attributes=modify_solver_attribute("verbose", false))
+ v = parsimonious_flux_balance_analysis_vec(model, Tulip.Optimizer; modifications=modify_constraint(glucose, -12, -12), qp_solver=modify_solver(OSQP.Optimizer), qp_solver_attributes=modify_solver_attribute("verbose", false))
+
+ @test isapprox(d["PGM"], -14.751259785795853, atol = 1e-3)
+ @test isapprox(v[8], -14.751259785795853, atol = 1e-3)
end
diff --git a/test/analysis/utilities.jl b/test/analysis/utilities.jl
new file mode 100644
index 0000000000000000000000000000000000000000..fe3ac19a3f158f3e7783305527d5224ef7fca7f2
--- /dev/null
+++ b/test/analysis/utilities.jl
@@ -0,0 +1,44 @@
+@testset "Analysis utilities" begin
+ model = read_model(
+ download_data_file(
+ "http://bigg.ucsd.edu/static/models/e_coli_core.json",
+ joinpath("data", "e_coli_core.json"),
+ "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
+ ),
+ )
+ @test length(model.reactions) == 95 # read in correctly
+
+ # FBA
+ biomass = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")
+
+ cons = Dict("EX_glc__D_e" => (-12.0, -12.0))
+ optimizer = Tulip.Optimizer # quiet by default
+ sol = flux_balance_analysis_dict(model, optimizer; modifications=modify_objective(biomass))
+
+ # atom tracker
+ atom_fluxes = atom_exchange(sol, model)
+ @test isapprox(atom_fluxes["C"], -37.1902, atol = 1e-3)
+
+ # exchange trackers
+ consuming, producing = exchange_reactions(sol; verbose = false)
+ @test isapprox(consuming["EX_nh4_e"], -4.76532, atol = 1e-3)
+
+ # metabolite trackers
+ consuming, producing = metabolite_fluxes(sol, model)
+ @test isapprox(consuming["atp_c"]["PFK"], -7.47738, atol = 1e-3)
+ @test isapprox(producing["atp_c"]["PYK"], 1.75818, atol = 1e-3)
+
+ # set bounds
+ cbm = make_optimization_model(model, optimizer)
+ ubs = cbm[:ubs]
+ lbs = cbm[:lbs]
+ glucose_index = model[findfirst(model.reactions, "EX_glc__D_e")]
+ o2_index = model[findfirst(model.reactions, "EX_o2_e")]
+ atpm_index = model[findfirst(model.reactions, "ATPM")]
+ set_bound(glucose_index, cbm; ub = -1.0, lb = -1.0)
+ @test normalized_rhs(ubs[glucose_index]) == -1.0
+ @test normalized_rhs(lbs[glucose_index]) == 1.0
+ set_bound(o2_index, cbm; ub = 1.0, lb = 1.0)
+ @test normalized_rhs(ubs[o2_index]) == 1.0
+ @test normalized_rhs(lbs[o2_index]) == -1.0
+end
diff --git a/test/base/solver.jl b/test/base/solver.jl
index 5f7f8ab7954e6657fb8cee911e84586d2d61dd30..bdfa647824eb087cdc2569542704fe9fba36fe26 100644
--- a/test/base/solver.jl
+++ b/test/base/solver.jl
@@ -3,15 +3,13 @@
cp = test_simpleLP()
optimizer = GLPK.Optimizer
lp = optimize_model(cp, optimizer)
- x = lp[:x]
@test termination_status(lp) === MOI.OPTIMAL
- sol = COBREXA.JuMP.value.(x)
+ sol = JuMP.value.(lp[:x])
@test sol ≈ [1.0, 2.0]
optimizer = Clp.Optimizer
lp = optimize_model(cp, optimizer)
- x = lp[:x]
@test termination_status(lp) === MOI.OPTIMAL
- sol = COBREXA.JuMP.value.(x)
+ sol = JuMP.value.(lp[:x])
@test sol ≈ [1.0, 2.0]
end
diff --git a/test/reconstruction/model_manipulations_test.jl b/test/reconstruction/model_manipulations_test.jl
index 8822ae7dc56e1d50c7bfffe6e8afc5c59503ce55..a3b74cdc33ccf22492e514a57b758ece67aca1b5 100644
--- a/test/reconstruction/model_manipulations_test.jl
+++ b/test/reconstruction/model_manipulations_test.jl
@@ -26,7 +26,7 @@
rxns = [r1, r2]
- model = CobraModel()
+ model = StandardModel()
model.id = "model"
model.reactions = rxns
model.metabolites = mets
diff --git a/test/sampling/sampling_tools_test.jl b/test/sampling/sampling_tools_test.jl
index 92611296bdeb91fd8fa33f25c5cbdc6e0eee107c..dea4e51d7c05aeab1f43c797fe2f0a3ade0ea62f 100644
--- a/test/sampling/sampling_tools_test.jl
+++ b/test/sampling/sampling_tools_test.jl
@@ -8,28 +8,21 @@
),
)
- optimizer = COBREXA.Tulip.Optimizer
biomass = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")
- cons = Dict("EX_glc__D_e" => (-12.0, -12.0))
- atts = Dict("IPM_IterationsLimit" => 110)
-
- sol = fba(model, optimizer; objective_func = biomass, constraints = cons)
-
- cons["BIOMASS_Ecoli_core_w_GAM"] =
- (sol["BIOMASS_Ecoli_core_w_GAM"] * 0.99, sol["BIOMASS_Ecoli_core_w_GAM"])
+ glucose = findfirst(model.reactions, "EX_glc__D_e")
+ cbm = flux_balance_analysis(model, Tulip.Optimizer; modifications=[modify_objective(biomass), modify_constraint(glucose, -12, -12), modify_solver_attribute("IPM_IterationsLimit", 500)])
+ biomass_index = model[biomass]
+ λ = JuMP.value(cbm[:x][biomass_index])
+ modify_constraint(biomass, 0.99*λ, λ)(model, cbm)
samples = hit_and_run(
100_000,
- model,
- optimizer;
+ cbm,
keepevery = 10,
samplesize = 5000,
- constraints = cons,
- solver_attributes = atts,
)
- # TODO
- # @test isapprox(mean(samples[64, :]), 8.9, atol = 0.5) # only tests if the sampler approximately converged
+ @test isapprox(mean(samples[64, :]), 8.9, atol = 0.5) # only tests if the sampler approximately converged
# samples = achr(
# 100_000,
diff --git a/test/types/model_test.jl b/test/types/model_test.jl
index 579077905bd2ee4af3dde541bc5a8fe111f02d1b..2c5b7958d3cf9a8543475419cde150d04224621d 100644
--- a/test/types/model_test.jl
+++ b/test/types/model_test.jl
@@ -31,7 +31,7 @@
genes = [g1, g2, g3]
rxns = [r1, r2, r3, r4]
- model = CobraModel()
+ model = StandardModel()
model.id = "model"
model.reactions = rxns
model.metabolites = mets