Commit 50735009 authored by Miroslav Kratochvil's avatar Miroslav Kratochvil
Browse files

remove unused code, sprinkle TODOs

parent a9677a09
Pipeline #55043 failed with stages
in 9 minutes and 32 seconds
......@@ -382,86 +382,3 @@ function Base.convert(::Type{StandardModel}, model::MetabolicModel)
genes = modelgenes,
)
end
#TODO generalize these to other model types
"""
reaction_bounds(model::StandardModel, rid::String)
Return lower and upper bounds for `rid` in `model`.
"""
function reaction_bounds(model::StandardModel, rid::String)
model.reactions[rid].lb, model.reactions[rid].ub
end
"""
is_reaction_reversible(model::StandardModel, rid::String)
Check if reaction `rid` in `model` is reversible.
"""
function is_reaction_reversible(model::StandardModel, rid::String)
lb, ub = reaction_bounds(model, rid)
lb < 0 && ub > 0
end
"""
is_reaction_forward_only(model::StandardModel, rid::String)
Check if reaction `rid` in `model` is forward only.
"""
function is_reaction_forward_only(model::StandardModel, rid::String)
lb, ub = reaction_bounds(model, rid)
lb >= 0 && ub > 0
end
"""
is_reaction_backward_only(model::StandardModel, rid::String)
Check if reaction `rid` in `model` is backward only.
"""
function is_reaction_backward_only(model::StandardModel, rid::String)
lb, ub = reaction_bounds(model, rid)
lb < 0 && ub <= 0
end
"""
is_reaction_unidirectional(model::StandardModel, rid::String)
Check if reaction `rid` in `model` is unidirectional.
"""
function is_reaction_unidirectional(model::StandardModel, rid::String)
is_reaction_forward_only(model, rid) || is_reaction_backward_only(model, rid)
end
"""
is_reaction_blocked(model::StandardModel, rid::String)
Check if reaction `rid` in `model` is blocked.
"""
function is_reaction_blocked(model::StandardModel, rid::String)
lb, ub = reaction_bounds(model, rid)
lb == ub == 0
end
"""
reaction_has_multiple_isozymes(model::StandardModel, rid::String)
Check if reaction `rid` in `model` is catalyzed by multiple enzymes,
i.e. it has isozymes according to the gene reaction rules.
"""
function reaction_has_multiple_isozymes(model::StandardModel, rid::String)
length(reaction_gene_association(model, rid)) > 1
end
"""
reaction_has_valid_gene_association(model::StandardModel, rid::String)
Check if reaction `rid` in `model` has a gene reaction rule entry.
"""
function reaction_has_valid_gene_association(model::StandardModel, rid::String)
#TODO simplify this once COBREXA enforces universal rules for GRR representation
haskey(model.reactions, rid) &&
!isnothing(reaction_gene_association(model, rid)) &&
reaction_gene_association(model, rid) != [[]] &&
!isempty(first(reaction_gene_association(model, rid)))
end
......@@ -39,3 +39,5 @@ protein_mass_group_dict(model::GeckoModel, opt_model) =
A pipe-able variant of [`mass_group_dict`](@ref).
"""
protein_mass_group_dict(model::GeckoModel) = x -> mass_group_dict(model, x)
#TODO: total protein mass for sMOMENT
"""
_build_irreversible_stoichiometric_matrix(model::StandardModel)
Return a stoichiometric matrix where all reactions are forward only i.e. only
positive fluxes are allowed. To accomplish this for models with isozymes,
so-called arm reactions are included. Note, reactions that are irreversible
in the original model will be irreversible in this model. E.g., if a reaction
is forward only in the original model, then there will be no reverse component
for this reaction in the irreversible stoichiometric matrix.
"""
function _build_irreversible_stoichiometric_matrix(
model::StandardModel,
rid_isozymes = Dict{String,Vector{Isozyme}}(),
)
# components used to build stoichiometric matrix
S_components = ( #TODO add size hints if possible
row_idxs = Vector{Int}(),
col_idxs = Vector{Int}(),
coeffs = Vector{Float64}(),
lbs = Vector{Float64}(),
ubs = Vector{Float64}(),
)
# establish the ordering in a named tuple
idxs = ( #: pseudo metabolites and reactions are added to model
met_idxs = Dict{String,Int}(),
rxn_idxs = Dict{String,Int}(),
max_rxn_idx = [1], #TODO maybe fix, this is a dodgy way of adding a counter to a named tuple
max_met_idx = [1], #TODO maybe fix, this is a dodgy way of adding a counter to a named tuple
pseudo_met_idx = [1], #TODO maybe fix, this is a dodgy way of adding a counter to a named tuple
)
#TODO for the counter thing, basically I wanted e.g. max_rxn_idx = 1 and then update it,
#TODO but named tuples are immutable... :(
# fill the matrix entries
#: blocked treated as reversible because unclear what direction the reaction would go
for rid in reactions(model)
if haskey(rid_isozymes, rid) && length(rid_isozymes[rid]) > 1
if is_reaction_unidirectional(model, rid)
dir = is_reaction_forward_only(model, rid) ? "§FOR" : "§REV"
_add_isozyme_to_irrev_stoich_mat(
model,
rid_isozymes,
rid,
idxs,
S_components,
dir,
)
elseif is_reaction_reversible(model, rid) || is_reaction_blocked(model, rid)
_add_isozyme_to_irrev_stoich_mat(
model,
rid_isozymes,
rid,
idxs,
S_components,
"§FOR",
)
_add_isozyme_to_irrev_stoich_mat(
model,
rid_isozymes,
rid,
idxs,
S_components,
"§REV",
)
else
@warn "Unhandled bound type for $rid"
end
else # no grr or single enzyme only
if is_reaction_unidirectional(model, rid)
dir = is_reaction_forward_only(model, rid) ? "§FOR" : "§REV"
_add_enzyme_to_irrev_stoich_mat(model, rid, idxs, S_components, dir)
elseif is_reaction_reversible(model, rid) || is_reaction_blocked(model, rid)
_add_enzyme_to_irrev_stoich_mat(model, rid, idxs, S_components, "§FOR")
_add_enzyme_to_irrev_stoich_mat(model, rid, idxs, S_components, "§REV")
else
@warn "Unhandled bound type for $rid"
end
end
end
S = sparse(
S_components.row_idxs,
S_components.col_idxs,
S_components.coeffs,
length(idxs.met_idxs),
length(idxs.rxn_idxs),
)
return S, S_components.lbs, S_components.ubs, idxs.rxn_idxs, idxs.met_idxs
end
"""
_add_enzyme_to_irrev_stoich_mat(model, rid, idxs, S_components, dir)
Add entries to the components that will be used to build the stoichiometric
matrix. Simple variant that does not deal with isozymes and arm reactions.
"""
function _add_enzyme_to_irrev_stoich_mat(model::StandardModel, rid, idxs, S_components, dir)
idxs.rxn_idxs[rid*dir] = idxs.max_rxn_idx[1]
idxs.max_rxn_idx[1] += 1
fix_sign = dir == "§FOR" ? 1 : -1 # change direction of reaction
for (mid, coeff) in reaction_stoichiometry(model, rid)
if !haskey(idxs.met_idxs, mid)
idxs.met_idxs[mid] = idxs.max_met_idx[1]
idxs.max_met_idx[1] += 1
end
push!(S_components.row_idxs, idxs.met_idxs[mid])
push!(S_components.col_idxs, idxs.rxn_idxs[rid*dir])
push!(S_components.coeffs, fix_sign * coeff)
end
lb, ub = abs.(reaction_bounds(model, rid)) # assumes lb < ub
if dir == "§FOR"
is_reaction_reversible(model, rid) ? push!(S_components.lbs, 0) :
push!(S_components.lbs, lb)
push!(S_components.ubs, ub)
else
is_reaction_reversible(model, rid) ? push!(S_components.lbs, 0) :
push!(S_components.lbs, ub)
push!(S_components.ubs, lb)
end
end
"""
_add_isozyme_to_irrev_stoich_mat(
model::StandardModel,
rid,
idxs,
S_components,
dir,
)
Add entries to the components that will be used to build the stoichiometric matrix.
Complex variant that deals with isozymes and arm reactions.
"""
function _add_isozyme_to_irrev_stoich_mat(
model::StandardModel,
rid_isoyzmes,
rid,
idxs,
S_components,
dir,
)
# add pseudo metabolite
pm = "§PM$(idxs.pseudo_met_idx[1])"
idxs.pseudo_met_idx[1] += 1
idxs.met_idxs[pm] = idxs.max_met_idx[1]
idxs.max_met_idx[1] += 1
# find half reactions to get arm reaction
lhs = []
rhs = []
for (mid, coeff) in reaction_stoichiometry(model, rid)
if !haskey(idxs.met_idxs, mid)
idxs.met_idxs[mid] = idxs.max_met_idx[1]
idxs.max_met_idx[1] += 1
end
if coeff <= 0
push!(lhs, (mid, coeff))
else
push!(rhs, (mid, coeff))
end
end
product_half_reaction = dir == "§FOR" ? rhs : lhs
reagent_half_reaction = dir == "§FOR" ? lhs : rhs
# add arm reaction
fix_sign = dir == "§FOR" ? 1 : -1 # change direction of reaction
pr = rid * "§ARM" * dir
idxs.rxn_idxs[pr] = idxs.max_rxn_idx[1] #! this needs to get added first because of blocked possibility
idxs.max_rxn_idx[1] += 1
push!(S_components.row_idxs, idxs.met_idxs[pm])
push!(S_components.col_idxs, idxs.rxn_idxs[pr])
push!(S_components.coeffs, 1)
for (mid, coeff) in reagent_half_reaction
push!(S_components.row_idxs, idxs.met_idxs[mid])
push!(S_components.col_idxs, idxs.rxn_idxs[pr])
push!(S_components.coeffs, fix_sign * coeff)
end
# add bounds for ARM reaction that corresponds to original model's bounds
lb, ub = abs.(reaction_bounds(model, rid)) # assumes lb < ub
if dir == "§FOR"
is_reaction_reversible(model, rid) ? push!(S_components.lbs, 0) :
push!(S_components.lbs, lb)
push!(S_components.ubs, ub)
else
is_reaction_reversible(model, rid) ? push!(S_components.lbs, 0) :
push!(S_components.lbs, ub)
push!(S_components.ubs, lb)
end
# add isozyme reactions
for (i, _) in enumerate(rid_isoyzmes[rid])
iso_rid = rid * "§ISO$i" * dir
idxs.rxn_idxs[iso_rid] = idxs.max_rxn_idx[1]
idxs.max_rxn_idx[1] += 1
push!(S_components.row_idxs, idxs.met_idxs[pm])
push!(S_components.col_idxs, idxs.rxn_idxs[iso_rid])
push!(S_components.coeffs, -1)
for (mid, coeff) in product_half_reaction
push!(S_components.row_idxs, idxs.met_idxs[mid])
push!(S_components.col_idxs, idxs.rxn_idxs[iso_rid])
push!(S_components.coeffs, fix_sign * coeff)
end
# add bounds
push!(S_components.lbs, 0)
if is_reaction_blocked(model, rid)
push!(S_components.ubs, 0)
else
push!(S_components.ubs, 1000) # arbitrary upper bound
end
end
end
"""
_order_id_to_idx_dict(id_to_idx_dict)
Return the keys of `id_to_idx_dict` sorted by the values, which
are taken to be the indices. This is a helper function for
[`reactions`](@ref) and [`metabolites`](@ref).
"""
function _order_id_to_idx_dict(dmap)
ks = collect(keys(dmap))
vs = collect(values(dmap))
return ks[sortperm(vs)]
end
......@@ -19,57 +19,7 @@ permanently irreversible in the model, i.e. changing their bounds to make them
reversible will have no effect.
"""
function change_bound!(model::GeckoModel, id; lower = nothing, upper = nothing)
gene_idx = first(indexin([id], model.gene_ids))
if isnothing(gene_idx)
flux_for_idx = findfirst(
x -> x == id * "§ARM§FOR" || x == id * "§FOR",
model.irrev_reaction_ids,
)
if !isnothing(flux_for_idx)
if !isnothing(lower)
if lower <= 0
model.xl[flux_for_idx] = 0
else
model.xl[flux_for_idx] = lower
end
end
if !isnothing(upper)
if upper <= 0
model.xu[flux_for_idx] = 0
else
model.xu[flux_for_idx] = upper
end
end
end
flux_rev_idx = findfirst(
x -> x == id * "§ARM§REV" || x == id * "§REV",
model.irrev_reaction_ids,
)
if !isnothing(flux_rev_idx)
if !isnothing(lower)
if lower >= 0
model.xu[flux_rev_idx] = 0
else
model.xu[flux_rev_idx] = -lower
end
if !isnothing(upper)
if upper >= 0
model.xl[flux_rev_idx] = 0
else
model.xl[flux_rev_idx] = -upper
end
end
end
end
else
n = length(model.irrev_reaction_ids)
!isnothing(lower) && (model.xl[n+gene_idx] = lower)
!isnothing(upper) && (model.xu[n+gene_idx] = upper)
end
return nothing
#TODO
end
"""
......@@ -84,7 +34,5 @@ function change_bounds!(
lower = fill(nothing, length(ids)),
upper = fill(nothing, length(ids)),
)
for (id, lower, upper) in zip(ids, lower, upper)
change_bound!(model, id; lower = lower, upper = upper)
end
#TODO
end
......@@ -19,47 +19,7 @@ permanently irreversible in the model, i.e. changing their bounds to make them
reversible will have no effect.
"""
function change_bound!(model::SMomentModel, id; lower = nothing, upper = nothing)
flux_for_idx =
findfirst(x -> x == id * "§ARM§FOR" || x == id * "§FOR", model.irrev_reaction_ids)
if !isnothing(flux_for_idx)
if !isnothing(lower)
if lower <= 0
model.xl[flux_for_idx] = 0
else
model.xl[flux_for_idx] = lower
end
end
if !isnothing(upper)
if upper <= 0
model.xu[flux_for_idx] = 0
else
model.xu[flux_for_idx] = upper
end
end
end
flux_rev_idx =
findfirst(x -> x == id * "§ARM§REV" || x == id * "§REV", model.irrev_reaction_ids)
if !isnothing(flux_rev_idx)
if !isnothing(lower)
if lower >= 0
model.xu[flux_rev_idx] = 0
else
model.xu[flux_rev_idx] = -lower
end
if !isnothing(upper)
if upper >= 0
model.xl[flux_rev_idx] = 0
else
model.xl[flux_rev_idx] = -upper
end
end
end
end
return nothing
# TODO
end
"""
......@@ -74,7 +34,5 @@ function change_bounds!(
lower = fill(nothing, length(ids)),
upper = fill(nothing, length(ids)),
)
for (id, lower, upper) in zip(ids, lower, upper)
change_bound!(model, id; lower = lower, upper = upper)
end
# TODO
end
......@@ -37,16 +37,16 @@ end
# set up the workers for Distributed, so that the tests that require more
# workers do not unnecessarily load the stuff multiple times
#W = addprocs(2)
#t = @elapsed @everywhere using COBREXA, Tulip, JuMP
#print_timing("import of packages", t)
#t = @elapsed @everywhere begin
#model = Model(Tulip.Optimizer)
#@variable(model, 0 <= x <= 1)
#@objective(model, Max, x)
#optimize!(model)
#end
#print_timing("JuMP+Tulip code warmup", t)
W = addprocs(2)
t = @elapsed @everywhere using COBREXA, Tulip, JuMP
print_timing("import of packages", t)
t = @elapsed @everywhere begin
model = Model(Tulip.Optimizer)
@variable(model, 0 <= x <= 1)
@objective(model, Max, x)
optimize!(model)
end
print_timing("JuMP+Tulip code warmup", t)
# make sure there's a directory for temporary data
tmpdir = "tmpfiles"
......@@ -59,17 +59,16 @@ run_test_file("data_downloaded.jl")
# import base files
@testset "COBREXA test suite" begin
#run_test_dir(joinpath("base", "types", "abstract"), "Abstract types")
#run_test_dir(joinpath("base", "types"), "Base model types")
#run_test_dir(joinpath("base", "logging"), "Logging")
#run_test_dir("base", "Base functionality")
#run_test_dir(joinpath("base", "utils"), "Utilities")
#run_test_dir("io", "I/O functions")
#run_test_dir("reconstruction")
#run_test_dir("analysis")
#run_test_dir(joinpath("analysis", "sampling"), "Sampling")
#run_test_file("analysis/smoment.jl")
run_test_file("analysis/gecko.jl")
run_test_dir(joinpath("base", "types", "abstract"), "Abstract types")
run_test_dir(joinpath("base", "types"), "Base model types")
run_test_dir(joinpath("base", "logging"), "Logging")
run_test_dir("base", "Base functionality")
run_test_dir(joinpath("base", "utils"), "Utilities")
run_test_dir("io", "I/O functions")
run_test_dir("reconstruction")
run_test_dir("analysis")
run_test_dir(joinpath("analysis", "sampling"), "Sampling")
run_test_file("aqua.jl")
end
#rmprocs(W)
rmprocs(W)
Markdown is supported
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