Unverified Commit 02ae1d62 authored by Miroslav Kratochvil's avatar Miroslav Kratochvil Committed by GitHub
Browse files

Merge pull request #515 from LCSB-BioCore/develop

develop→master for 1.1 release
parents 37dab779 0cb74ae7
Pipeline #49629 passed with stages
in 21 minutes and 33 seconds
name = "COBREXA"
uuid = "babc4406-5200-4a30-9033-bf5ae714c842"
authors = ["The developers of COBREXA.jl"]
version = "1.0.6"
version = "1.1"
[deps]
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
......@@ -40,12 +40,12 @@ julia = "1"
Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
OSQP = "ab2f91bb-94b4-55e3-9ba0-7f65df51de79"
SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"
GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
[targets]
test = ["Aqua", "Downloads", "OSQP", "SHA", "Test", "Tulip", "GLPK"]
......@@ -73,8 +73,8 @@ function COBREXA.objective(m::CircularModel)
end
COBREXA.bounds(m::CircularModel) = (
spzeros(n_reactions(m)), # lower bounds
sparse(ones(n_reactions(m))), # upper bounds
zeros(n_reactions(m)), # lower bounds
ones(n_reactions(m)), # upper bounds
)
function COBREXA.stoichiometry(m::CircularModel)
......
......@@ -5,28 +5,28 @@ A variant of FBA that returns a vector of fluxes in the same order as reactions
of the model, if the solution is found.
Arguments are passed to [`flux_balance_analysis`](@ref).
This function is kept for backwards compatibility, use [`flux_vector`](@ref)
instead.
"""
function flux_balance_analysis_vec(args...; kwargs...)::Union{Vector{Float64},Nothing}
optmodel = flux_balance_analysis(args...; kwargs...)
is_solved(optmodel) || return nothing
value.(optmodel[:x])
end
flux_balance_analysis_vec(args...; kwargs...)::Maybe{Vector{Float64}} =
flux_vector(flux_balance_analysis(args...; kwargs...))
"""
flux_balance_analysis_dict(model::MetabolicModel, args...)::Union{Dict{String, Float64},Nothing}
A variant of FBA that returns a dictionary assigning fluxes to reactions, if
the solution is found. Arguments are passed to [`flux_balance_analysis`](@ref).
This function is kept for backwards compatibility, use [`flux_dict`](@ref)
instead.
"""
function flux_balance_analysis_dict(
flux_balance_analysis_dict(
model::MetabolicModel,
args...;
kwargs...,
)::Union{Dict{String,Float64},Nothing}
v = flux_balance_analysis_vec(model, args...; kwargs...)
isnothing(v) && return nothing
Dict(zip(reactions(model), v))
end
)::Maybe{Dict{String,Float64}} =
flux_dict(model, flux_balance_analysis(model, args...; kwargs...))
"""
flux_balance_analysis(
......
......@@ -140,18 +140,11 @@ mins, maxs = flux_variability_analysis_dict(
```
"""
function flux_variability_analysis_dict(model::MetabolicModel, optimizer; kwargs...)
vs = flux_variability_analysis(
model,
optimizer;
kwargs...,
ret = m -> JuMP.value.(m[:x]),
)
fluxes = flux_variability_analysis(model, optimizer; kwargs..., ret = flux_vector)
rxns = reactions(model)
dicts = zip.(Ref(rxns), fluxes)
return (
Dict(zip(rxns, [Dict(zip(rxns, fluxes)) for fluxes in vs[:, 1]])),
Dict(zip(rxns, [Dict(zip(rxns, fluxes)) for fluxes in vs[:, 2]])),
)
return (Dict(rxns .=> Dict.(dicts[:, 1])), Dict(rxns .=> Dict.(dicts[:, 2])))
end
"""
......@@ -166,9 +159,5 @@ function _max_variability_flux(opt_model, rid, ret)
@objective(opt_model, sense, var)
optimize!(opt_model)
if is_solved(opt_model)
return ret(opt_model)
else
return nothing
end
is_solved(opt_model) ? ret(opt_model) : nothing
end
......@@ -69,7 +69,7 @@ minimize_metabolic_adjustment(flux_ref::Vector{Float64}) =
end
"""
minimize_metabolic_adjustment(flux_ref_dict::Dict{String, Float64})
minimize_metabolic_adjustment(flux_ref_dict::Dict{String, Float64})
Overload of [`minimize_metabolic_adjustment`](@ref) that works with a
dictionary of fluxes.
......@@ -87,30 +87,22 @@ minimize_metabolic_adjustment(flux_ref_dict::Dict{String,Float64}) =
Perform minimization of metabolic adjustment (MOMA) and return a vector of fluxes in the
same order as the reactions in `model`. Arguments are forwarded to
[`minimize_metabolic_adjustment`](@ref) internally.
"""
function minimize_metabolic_adjustment_analysis_vec(args...; kwargs...)
opt_model = minimize_metabolic_adjustment_analysis(args...; kwargs...)
isnothing(opt_model) && return nothing
return value.(opt_model[:x])
end
This function is kept for backwards compatibility, use [`flux_vector`](@ref)
instead.
"""
minimize_metabolic_adjustment_analysis_vec(args...; kwargs...) =
flux_vector(minimize_metabolic_adjustment_analysis(args...; kwargs...))
"""
minimize_metabolic_adjustment_analysis_dict(args...; kwargs...)
minimize_metabolic_adjustment_analysis_dict(model::MetabolicModel, args...; kwargs...)
Perform minimization of metabolic adjustment (MOMA) and return a dictionary mapping the
reaction IDs to fluxes. Arguments are forwarded to [`minimize_metabolic_adjustment`](@ref)
internally.
"""
function minimize_metabolic_adjustment_analysis_dict(
model::MetabolicModel,
args...;
kwargs...,
)
opt_fluxes = minimize_metabolic_adjustment_analysis_vec(model, args...; kwargs...)
isnothing(opt_fluxes) && return nothing
return Dict(zip(reactions(model), opt_fluxes))
end
This function is kept for backwards compatibility, use [`flux_vector`](@ref)
instead.
"""
minimize_metabolic_adjustment_analysis_dict(model::MetabolicModel, args...; kwargs...) =
flux_dict(model, minimize_metabolic_adjustment_analysis(model, args...; kwargs...))
......@@ -29,18 +29,10 @@ that active site number and unit issues are prevented.
# Example
```
sol = flux_balance_analysis_dict(
model,
Tulip.Optimizer;
modifications = [
add_moment_constraints(
ksas,
protein_mass_fraction
),
change_constraint("EX_glc__D_e", lb = -1000),
],
flux_balance_analysis(
...,
modifications = [ add_moment_constraints(my_kcats, 0.6) ],
)
```
"""
add_moment_constraints(kcats::Dict{String,Float64}, protein_mass_fraction::Float64) =
(model, opt_model) -> begin
......
......@@ -100,14 +100,12 @@ Perform parsimonious flux balance analysis on `model` using `optimizer`.
Returns a vector of fluxes in the same order as the reactions in `model`.
Arguments are forwarded to [`parsimonious_flux_balance_analysis`](@ref)
internally.
"""
function parsimonious_flux_balance_analysis_vec(args...; kwargs...)
opt_model = parsimonious_flux_balance_analysis(args...; kwargs...)
isnothing(opt_model) && return nothing
return value.(opt_model[:x])
end
This function is kept for backwards compatibility, use [`flux_vector`](@ref)
instead.
"""
parsimonious_flux_balance_analysis_vec(args...; kwargs...) =
flux_vector(parsimonious_flux_balance_analysis(args...; kwargs...))
"""
parsimonious_flux_balance_analysis_dict(model::MetabolicModel, args...; kwargs...)
......@@ -115,11 +113,9 @@ end
Perform parsimonious flux balance analysis on `model` using `optimizer`.
Returns a dictionary mapping the reaction IDs to fluxes. Arguments are
forwarded to [`parsimonious_flux_balance_analysis`](@ref) internally.
"""
function parsimonious_flux_balance_analysis_dict(model::MetabolicModel, args...; kwargs...)
opt_fluxes = parsimonious_flux_balance_analysis_vec(model, args...; kwargs...)
isnothing(opt_fluxes) && return nothing
return Dict(zip(reactions(model), opt_fluxes))
end
This function is kept for backwards compatibility, use [`flux_dict`](@ref)
instead.
"""
parsimonious_flux_balance_analysis_dict(model::MetabolicModel, args...; kwargs...) =
flux_dict(model, parsimonious_flux_balance_analysis(model, args...; kwargs...))
......@@ -32,29 +32,23 @@ function make_optimization_model(model::MetabolicModel, optimizer; sense = MOI.M
end
"""
is_solved(optmodel)
is_solved(opt_model)
Return `true` if `optmodel` solved successfully (solution is optimal or locally
Return `true` if `opt_model` solved successfully (solution is optimal or locally
optimal). Return `false` if any other termination status is reached.
Termination status is defined in the documentation of `JuMP`.
"""
function is_solved(optmodel)
termination_status(optmodel) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ? true : false
end
is_solved(opt_model) = termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
"""
optimize_objective(optmodel)::Union{Float64,Nothing}
optimize_objective(opt_model)::Maybe{Float64}
Shortcut for running JuMP `optimize!` on a model and returning the objective
value, if solved.
"""
function optimize_objective(optmodel)::Maybe{Float64}
optimize!(optmodel)
if is_solved(optmodel)
objective_value(optmodel)
else
nothing
end
function optimize_objective(opt_model)::Maybe{Float64}
optimize!(opt_model)
solved_objective_value(opt_model)
end
"""
......@@ -88,3 +82,43 @@ function set_optmodel_bound!(
isnothing(lb) || set_normalized_rhs(opt_model[:lbs][vidx], -lb)
isnothing(ub) || set_normalized_rhs(opt_model[:ubs][vidx], ub)
end
"""
solved_objective_value(opt_model)::Maybe{Float64}
Returns the current objective value of a model, if solved.
# Example
```
solved_objective_value(flux_balance_analysis(model, ...))
```
"""
solved_objective_value(opt_model)::Maybe{Float64} =
is_solved(opt_model) ? objective_value(opt_model) : nothing
"""
flux_vector(opt_model)::Maybe{Vector{Float64}}
Returns a vector of fluxes of the model, if solved.
# Example
```
flux_vector(flux_balance_analysis(model, ...))
```
"""
flux_vector(opt_model)::Maybe{Vector{Float64}} =
is_solved(opt_model) ? value.(opt_model[:x]) : nothing
"""
flux_dict(model::MetabolicModel, opt_model)::Maybe{Dict{String, Float64}, Nothing}
Returns the fluxes of the model as a reaction-keyed dictionary, if solved.
# Example
```
flux_dict(model, flux_balance_analysis(model, ...))
```
"""
flux_dict(model::MetabolicModel, opt_model)::Maybe{Dict{String,Float64}} =
is_solved(opt_model) ? Dict(reactions(model) .=> value.(opt_model[:x])) : nothing
......@@ -14,28 +14,27 @@ mutable struct CoreModel <: MetabolicModel
S::SparseMat
b::SparseVec
c::SparseVec
xl::SparseVec
xu::SparseVec
xl::Vector{Float64}
xu::Vector{Float64}
rxns::Vector{String}
mets::Vector{String}
function CoreModel(
S::M,
b::V,
c::V,
xl::V,
xu::V,
rxns::K,
mets::K,
) where {V<:VecType,M<:MatType,K<:StringVecType}
S::MatType,
b::VecType,
c::VecType,
xl::VecType,
xu::VecType,
rxns::StringVecType,
mets::StringVecType,
)
all([length(b), length(mets)] .== size(S, 1)) ||
throw(DimensionMismatch("inconsistent number of metabolites"))
all([length(c), length(xl), length(xu), length(rxns)] .== size(S, 2)) ||
throw(DimensionMismatch("inconsistent number of reactions"))
new(sparse(S), sparse(b), sparse(c), sparse(xl), sparse(xu), rxns, mets)
new(sparse(S), sparse(b), sparse(c), collect(xl), collect(xu), rxns, mets)
end
end
......@@ -61,11 +60,11 @@ metabolites(a::CoreModel)::Vector{String} = a.mets
stoichiometry(a::CoreModel)::SparseMat = a.S
"""
bounds(a::CoreModel)::Tuple{SparseVec,SparseVec}
bounds(a::CoreModel)::Tuple{Vector{Float64},Vector{Float64}}
`CoreModel` flux bounds.
"""
bounds(a::CoreModel)::Tuple{SparseVec,SparseVec} = (a.xl, a.xu)
bounds(a::CoreModel)::Tuple{Vector{Float64},Vector{Float64}} = (a.xl, a.xu)
"""
balance(a::CoreModel)::SparseVec
......
......@@ -10,22 +10,16 @@ The linear model with additional coupling constraints in the form
mutable struct CoreModelCoupled <: MetabolicModel
lm::CoreModel
C::SparseMat
cl::SparseVec
cu::SparseVec
function CoreModelCoupled(
lm::MetabolicModel,
C::M,
cl::V,
cu::V,
) where {V<:VecType,M<:MatType}
cl::Vector{Float64}
cu::Vector{Float64}
function CoreModelCoupled(lm::MetabolicModel, C::MatType, cl::VecType, cu::VecType)
length(cu) == length(cl) ||
throw(DimensionMismatch("`cl` and `cu` need to have the same size"))
size(C) == (length(cu), n_reactions(lm)) ||
throw(DimensionMismatch("wrong dimensions of `C`"))
new(convert(CoreModel, lm), sparse(C), sparse(cl), sparse(cu))
new(convert(CoreModel, lm), sparse(C), collect(cl), collect(cu))
end
end
......@@ -92,11 +86,11 @@ The number of coupling constraints in a `CoreModelCoupled`.
n_coupling_constraints(a::CoreModelCoupled)::Int = size(a.C, 1)
"""
coupling_bounds(a::CoreModelCoupled)::Tuple{SparseVec,SparseVec}
coupling_bounds(a::CoreModelCoupled)::Tuple{Vector{Float64},Vector{Float64}}
Coupling bounds for a `CoreModelCoupled`.
"""
coupling_bounds(a::CoreModelCoupled)::Tuple{SparseVec,SparseVec} = (a.cl, a.cu)
coupling_bounds(a::CoreModelCoupled)::Tuple{Vector{Float64},Vector{Float64}} = (a.cl, a.cu)
"""
reaction_stoichiometry(model::CoreModelCoupled, rid::String)::Dict{String, Float64}
......
......@@ -47,7 +47,7 @@ larger than `large_flux_bound` are only stored if `keep_unbounded` is `true`.
# Example
```
julia> sol = flux_balance_analysis_dict(model, Tulip.Optimizer)
julia> sol = flux_dict(flux_balance_analysis(model, Tulip.Optimizer))
julia> fr = flux_summary(sol)
Biomass:
BIOMASS_Ecoli_core_w_GAM: 0.8739
......
......@@ -146,12 +146,8 @@ Get the bounds for reactions, assuming the information is stored in
`.lower_bound` and `.upper_bound`.
"""
bounds(model::JSONModel) = (
sparse([
get(rxn, "lower_bound", -_constants.default_reaction_bound) for rxn in model.rxns
]),
sparse([
get(rxn, "upper_bound", _constants.default_reaction_bound) for rxn in model.rxns
]),
[get(rxn, "lower_bound", -_constants.default_reaction_bound) for rxn in model.rxns],
[get(rxn, "upper_bound", _constants.default_reaction_bound) for rxn in model.rxns],
)
"""
......
......@@ -59,8 +59,8 @@ stoichiometry(m::MATModel) = sparse(m.mat["S"])
Extracts bounds from the MAT file, saved under `lb` and `ub`.
"""
bounds(m::MATModel) = (
sparse(reshape(get(m.mat, "lb", fill(-Inf, n_reactions(m), 1)), n_reactions(m))),
sparse(reshape(get(m.mat, "ub", fill(Inf, n_reactions(m), 1)), n_reactions(m))),
reshape(get(m.mat, "lb", fill(-Inf, n_reactions(m), 1)), n_reactions(m)),
reshape(get(m.mat, "ub", fill(Inf, n_reactions(m), 1)), n_reactions(m)),
)
"""
......
......@@ -47,12 +47,12 @@ function stoichiometry(model::SBMLModel)::SparseMat
end
"""
bounds(model::SBMLModel)::Tuple{SparseVec,SparseVec}
bounds(model::SBMLModel)::Tuple{Vector{Float64},Vector{Float64}}
Get the lower and upper flux bounds of model [`SBMLModel`](@ref). Throws `DomainError` in
case if the SBML contains mismatching units.
"""
function bounds(model::SBMLModel)::Tuple{SparseVec,SparseVec}
function bounds(model::SBMLModel)::Tuple{Vector{Float64},Vector{Float64}}
lbu, ubu = SBML.flux_bounds(model.sbml)
unit = lbu[1][2]
......@@ -67,7 +67,7 @@ function bounds(model::SBMLModel)::Tuple{SparseVec,SparseVec}
),
)
return sparse.((getvalue.(lbu), getvalue.(ubu)))
return (getvalue.(lbu), getvalue.(ubu))
end
"""
......
......@@ -144,29 +144,29 @@ function stoichiometry(model::StandardModel)::SparseMat
end
"""
lower_bounds(model::StandardModel)
lower_bounds(model::StandardModel)::Vector{Float64}
Return the lower bounds for all reactions in `model` in sparse format.
"""
lower_bounds(model::StandardModel)::SparseVec =
lower_bounds(model::StandardModel)::Vector{Float64} =
sparse([model.reactions[rxn].lb for rxn in reactions(model)])
"""
upper_bounds(model::StandardModel)
upper_bounds(model::StandardModel)::Vector{Float64}
Return the upper bounds for all reactions in `model` in sparse format.
Order matches that of the reaction ids returned in `reactions()`.
"""
upper_bounds(model::StandardModel)::SparseVec =
upper_bounds(model::StandardModel)::Vector{Float64} =
sparse([model.reactions[rxn].ub for rxn in reactions(model)])
"""
bounds(model::StandardModel)
bounds(model::StandardModel)::Tuple{Vector{Float64},Vector{Float64}}
Return the lower and upper bounds, respectively, for reactions in `model`.
Order matches that of the reaction ids returned in `reactions()`.
"""
bounds(model::StandardModel)::Tuple{SparseVec,SparseVec} =
bounds(model::StandardModel)::Tuple{Vector{Float64},Vector{Float64}} =
(lower_bounds(model), upper_bounds(model))
"""
......
......@@ -105,11 +105,11 @@ function stoichiometry(a::MetabolicModel)::SparseMat
end
"""
bounds(a::MetabolicModel)::Tuple{SparseVec,SparseVec}
bounds(a::MetabolicModel)::Tuple{Vector{Float64},Vector{Float64}}
Get the lower and upper flux bounds of a model.
"""
function bounds(a::MetabolicModel)::Tuple{SparseVec,SparseVec}
function bounds(a::MetabolicModel)::Tuple{Vector{Float64},Vector{Float64}}
_missing_impl_error(bounds, (a,))
end
......@@ -151,12 +151,12 @@ function n_coupling_constraints(a::MetabolicModel)::Int
end
"""
coupling_bounds(a::MetabolicModel)::Tuple{SparseVec,SparseVec}
coupling_bounds(a::MetabolicModel)::Tuple{Vector{Float64},Vector{Float64}}
Get the lower and upper bounds for each coupling bound in a model, as specified
by `coupling`. By default, the model does not have any coupling bounds.
"""
function coupling_bounds(a::MetabolicModel)::Tuple{SparseVec,SparseVec}
function coupling_bounds(a::MetabolicModel)::Tuple{Vector{Float64},Vector{Float64}}
return (spzeros(0), spzeros(0))
end
......
......@@ -26,8 +26,8 @@ function add_reactions!(model::CoreModel, rxns::Vector{Reaction})
Sadd = sparse(I, J, V, n_metabolites(model), length(rxns))
model.S = [model.S Sadd]
model.c = dropzeros([model.c; cs])
model.xu = sparse(ubs)
model.xl = sparse(lbs)
model.xu = ubs
model.xl = lbs
return nothing
end
......@@ -42,25 +42,25 @@ add_reaction!(model::CoreModel, rxn::Reaction) = add_reactions!(model, [rxn])
"""
add_reactions(
m::CoreModel,
s::V1,
b::V2,
s::VecType,
b::VecType,
c::AbstractFloat,
xl::AbstractFloat,
xu::AbstractFloat;
check_consistency = false,
) where {V1<:VecType,V2<:VecType}
)
Add reaction(s) to a `CoreModel` model `m`.
"""
function add_reactions(
m::CoreModel,
s::V1,
b::V2,
s::VecType,
b::VecType,
c::AbstractFloat,
xl::AbstractFloat,
xu::AbstractFloat;
check_consistency = false,
) where {V1<:VecType,V2<:VecType}
)
return add_reactions(
m,
sparse(reshape(s, (length(s), 1))),
......@@ -75,27 +75,27 @@ end
"""
add_reactions(
m::CoreModel,
s::V1,
b::V2,
s::VecType,
b::VecType,
c::AbstractFloat,
xl::AbstractFloat,
xu::AbstractFloat,
rxn::String,
mets::K;
check_consistency = false,
) where {V1<:VecType,V2<:VecType,K<:StringVecType}
)
"""
function add_reactions(
m::CoreModel,
s::V1,
b::V2,
s::VecType,
b::VecType,
c::AbstractFloat,
xl::AbstractFloat,
xu::AbstractFloat,
rxn::String,
mets::K;
mets::StringVecType;
check_consistency = false,
) where {V1<:VecType,V2<:VecType,K<:StringVecType}
)
return add_reactions(
m,
sparse(reshape(s, (length(s), 1))),
......@@ -112,23 +112,23 @@ end
"""
add_reactions(
m::CoreModel,
Sp::M,
b::V,
c::V,