Merge pull request #406 from LCSB-BioCore/mk-envelopes

implement objective envelopes (aka production envelopes)

src/analysis/envelopes.jl

+
+"""
+    envelope_lattice(model::MetabolicModel, rids::Vector{String}; kwargs...)
+
+Version of [`envelope_lattice`](@ref) that works on string reaction IDs instead
+of integer indexes.
+"""
+envelope_lattice(model::MetabolicModel, rids::Vector{String}; kwargs...) =
+    envelope_lattice(model, Vector{Int}(indexin(rids, reactions(model))); kwargs...)
+
+"""
+    envelope_lattice(
+        model::MetabolicModel,
+        ridxs::Vector{Int};
+        samples = 10,
+        ranges = collect(zip(bounds(model)...))[ridxs],
+        reaction_samples = fill(samples, length(ridxs)),
+    )
+
+Create a lattice (list of "tick" vectors) for reactions at indexes `ridxs` in a
+model. Arguments `samples`, `ranges`, and `reaction_samples` may be optionally
+specified to customize the latice creation process.
+"""
+envelope_lattice(
+    model::MetabolicModel,
+    ridxs::Vector{Int};
+    samples = 10,
+    ranges = collect(zip(bounds(model)...))[ridxs],
+    reaction_samples = fill(samples, length(ridxs)),
+) = (
+    lb .+ (ub - lb) .* ((1:s) .- 1) ./ max(s - 1, 1) for
+    (s, (lb, ub)) in zip(reaction_samples, ranges)
+)
+
+"""
+    objective_envelope(model::MetabolicModel, rids::Vector{String}, args...; kwargs...)
+
+Versioin of [`objective_envelope`](@ref) that works on string reaction IDs
+instead of integer indexes.
+"""
+objective_envelope(model::MetabolicModel, rids::Vector{String}, args...; kwargs...) =
+    objective_envelope(
+        model,
+        Vector{Int}(indexin(rids, reactions(model))),
+        args...;
+        kwargs...,
+    )
+
+"""
+    objective_envelope(
+        model::MetabolicModel,
+        ridxs::Vector{Int},
+        optimizer;
+        lattice = envelope_lattice(model, ridxs),
+        kwargs...,
+    )
+
+Compute an array of objective values for the `model` for rates of reactions
+specified `ridxs` fixed to a regular range of values between their respective
+lower and upper bounds.
+
+This can be used to compute a "production envelope" of a metabolite; but
+generalizes to any specifiable objective and to multiple dimensions of the
+examined space. To retrieve a production envelope of any metabolite, set the
+objective coefficient vector of the `model` to a vector that contains a single
+`1` for the exchange reaction that "outputs" this metabolite, and run
+[`objective_envelope`](@ref) with the exchange reaction of the "parameter"
+metabolite specified in `ridxs`.
+
+Returns a named tuple that contains `lattice` with reference values of the
+metabolites, and an N-dimensional array `values` with the computed objective
+values, where N is the number of specified reactions.  Because of the
+increasing dimensionality, the computation gets very voluminous with increasing
+length of `ridxs`.
+
+`kwargs` are internally forwarded to [`screen_optmodel_modifications`](@ref).
+"""
+objective_envelope(
+    model::MetabolicModel,
+    ridxs::Vector{Int},
+    optimizer;
+    lattice = envelope_lattice(model, ridxs),
+    kwargs...,
+) = (
+    lattice = collect.(lattice),
+    values = screen_optmodel_modifications(
+        model,
+        optimizer,
+        modifications = collect(
+            [(_, optmodel) -> begin
+                    for (i, ridx) in enumerate(ridxs)
+                        set_normalized_rhs(optmodel[:lbs][ridx], bounds[i])
+                        set_normalized_rhs(optmodel[:ubs][ridx], bounds[i])
+                    end
+                end] for bounds in Iterators.product(lattice...)
+        ),
+        analysis = screen_optimize_objective;
+        kwargs...,
+    ),
+)

src/analysis/screening.jl

@@ -143,3 +143,82 @@ A shortcut for [`screen`](@ref) that only works with model variants.
 """
 screen_variants(model, variants, analysis; workers = [myid()]) =
     screen(model; variants = variants, analysis = analysis, workers = workers)
+
+"""
+    screen_optimize_objective(_, optmodel)::Union{Float64,Nothing}
+
+A variant of [`optimize_objective`](@ref) directly usable in
+[`screen_optmodel_modifications`](@ref).
+"""
+screen_optimize_objective(_, optmodel)::Union{Float64,Nothing} =
+    optimize_objective(optmodel)
+
+"""
+    screen_optmodel_modifications(
+        model::MetabolicModel,
+        optimizer;
+        modifications::Array{V,N},
+        analysis = screen_optimize_objective,
+        workers = [myid()],
+    ) where {V<:AbstractVector,N}
+
+Screen multiple modifications of the same optimization model.
+
+This function is potentially more efficient than [`screen`](@ref) because it
+avoids making variants of the model structure and remaking of the optimization
+model. On the other hand, modification functions need to keep the optimization
+model in a recoverable state (one that leaves the model usable for the next
+modification), which limits the possible spectrum of modifications applied.
+
+Internally, `model` is distributed to `workers` and transformed into the
+optimization model using [`make_optimization_model`](@ref). With that,
+vectors of functions in `modifications` are consecutively applied, and the
+result of `analysis` function called on model are collected to an array of the
+same extent as `modifications`.
+
+Both the modification functions (in vectors) and the analysis function here
+have 2 parameters (as opposed to 1 with [`screen`](@ref)): first is the `model`
+(carried through as-is), second is the prepared JuMP optimization model that
+may be modified and acted upon. As an example, you can use modification
+[`change_constraint`](@ref) and analysis [`screen_optimize_objective`](@ref).
+"""
+function screen_optmodel_modifications(
+    model::MetabolicModel,
+    optimizer;
+    modifications::Array{V,N},
+    analysis = screen_optimize_objective,
+    workers = [myid()],
+) where {V<:AbstractVector,N}
+    save_model = :(
+        begin
+            local model = $model
+            $COBREXA.precache!(model)
+            (model, $COBREXA.make_optimization_model(model, $optimizer))
+        end
+    )
+    map(
+        fetch,
+        save_at.(workers, :cobrexa_screen_optmodel_modifications_data, Ref(save_model)),
+    )
+    save_model = nothing
+    map(fetch, save_at.(workers, :cobrexa_screen_optmodel_modifications_fn, Ref(analysis)))
+
+    res = dpmap(
+        mods -> :(
+            begin
+                local (model, optmodel) = cobrexa_screen_optmodel_modifications_data
+                for mod in $mods
+                    mod(model, optmodel)
+                end
+                cobrexa_screen_optmodel_modifications_fn(model, optmodel)
+            end
+        ),
+        CachingPool(workers),
+        modifications,
+    )
+
+    map(fetch, remove_from.(workers, :cobrexa_screen_optmodel_modifications_data))
+    map(fetch, remove_from.(workers, :cobrexa_screen_optmodel_modifications_fn))
+
+    return res
+end

src/base/solver.jl

@@ -59,6 +59,15 @@ function is_solved(optmodel)
     false
 end
 
+function optimize_objective(optmodel)::Union{Float64,Nothing}
+    optimize!(optmodel)
+    if is_solved(optmodel)
+        objective_value(optmodel)
+    else
+        nothing
+    end
+end
+
 """
     get_optmodel_bounds(opt_model)
 

test/analysis/envelopes.jl

+
+@testset "Envelopes" begin
+    m = load_model(model_paths["e_coli_core.xml"])
+
+    rxns = [1, 2, 3]
+
+    lat = collect.(envelope_lattice(m, rxns; samples = 3))
+    @test lat == [[0, 500, 1000], [-1000, 0, 1000], [-1000, 0, 1000]]
+    @test lat == collect.(envelope_lattice(m, reactions(m)[rxns]; samples = 3))
+
+    vals =
+        objective_envelope(
+            m,
+            reactions(m)[rxns],
+            Tulip.Optimizer;
+            lattice = lat,
+            workers = W,
+        ).values
+
+    @test size(vals) == (3, 3, 3)
+    @test count(isnothing, vals) == 15
+    @test isapprox(
+        filter(!isnothing, vals),
+        [
+            0.0,
+            0.0,
+            0.0,
+            0.704036947,
+            10.053282384,
+            10.053282365,
+            0.0,
+            0.0,
+            0.0,
+            0.873921506,
+            18.664485767,
+            20.412058183,
+        ],
+        atol = TEST_TOLERANCE,
+    )
+end