FVA warmup

src/analysis/sampling/warmup.jl

@@ -7,28 +7,53 @@
         workers = [myid()],
     )
 
-Generates warmup points for samplers by minimizing and maximizing random
-reactions. Very similar to
-[`flux_variability_analysis`](@ref).
+Generates FVA-like warmup points for samplers, by selecting random points by
+minimizing and maximizing reactions. Can not return more than 2 times the
+number of reactions in the model.
 """
-function warmup_from_variability(
-    n_points::Int,
-    model::MetabolicModel,
-    optimizer;
-    kwargs...
-)
+function warmup_from_variability(n_points::Int, model::MetabolicModel, optimizer; kwargs...)
+    nr = n_reactions(model)
+
+    n_points > 2 * nr && throw(
+        DomainError(
+            n_points,
+            "Variability method can not generate more than $(2*nr) points from this model",
+        ),
+    )
 
-    
+    sample = shuffle(vcat(1:nr, -(1:nr)))[begin:n_points]
+    warmup_from_variability(
+        -filter(x -> x < 0, sample),
+        filter(x -> x > 0, sample),
+        model,
+        optimizer;
+        kwargs...,
+    )
 end
 
+"""
+    warmup_from_variability(
+        min_reactions::Vector{Int},
+        max_reactions::Vector{Int},
+        model::MetabolicModel,
+        optimizer;
+        modifications=[],
+        workers::Vector{Int}=[myid()]
+    )::Matrix{Float64}
+
+Generate FVA-like warmup points for samplers, by minimizing and maximizing the
+specified reactions. The result is returned as a matrix, each point occupies as
+single column in the result.
+"""
 function warmup_from_variability(
     min_reactions::Vector{Int},
     max_reactions::Vector{Int},
     model::MetabolicModel,
     optimizer;
-    modifications=[],
-    workers::Vector{Int}=[myid()]
-)::Tuple{Matrix{Float64}, AbstractVector{Float64}, AbstractVector{Float64}}
+    modifications = [],
+    workers::Vector{Int} = [myid()],
+)::Matrix{Float64}
+
     # create optimization problem at workers, apply modifications
     save_model = :(
         begin
@@ -43,18 +68,20 @@ function warmup_from_variability(
 
     map(fetch, save_at.(workers, :cobrexa_sampling_warmup_optmodel, Ref(save_model)))
 
-    fluxes = dpmap(
-        rid -> :($COBREXA._FVA_optimize_reaction(
-            cobrexa_sampling_warmup_optmodel,
-            $rid,
-            optmodel -> value.(optmodel[:x]),
-        )),
-        CachingPool(workers),
-        vcat(-min_reactions, max_reactions),
+    fluxes = hcat(
+        dpmap(
+            rid -> :($COBREXA._FVA_optimize_reaction(
+                cobrexa_sampling_warmup_optmodel,
+                $rid,
+                optmodel -> $COBREXA.JuMP.value.(optmodel[:x]),
+            )),
+            CachingPool(workers),
+            vcat(-min_reactions, max_reactions),
+        )...,
     )
 
     # free the data on workers
     map(fetch, remove_from.(workers, :cobrexa_sampling_warmup_optmodel))
 
-    return hcat(fluxes...)
+    return fluxes
 end

test/analysis/sampling/warmup.jl

@@ -6,33 +6,16 @@
     )
     model = load_model(StandardModel, model_path)
 
-    # Serial test
-    ws, lbs, ubs = warmup(
+    pts = warmup_from_variability(
+        100,
         model,
         Tulip.Optimizer;
-        modifications = [change_constraint("EX_glc__D_e", -4, 4)],
-        warmup_points = collect(1:n_reactions(model)),
-        workerids = W,
+        modifications = [change_constraint("EX_glc__D_e", -2, 2)],
+        workers = W,
     )
 
-    ind = first(indexin(["EX_glc__D_e"], reactions(model)))
-    @test size(ws) == (95, 2)
-    @test size(ws[1, 1]) == (95,)
-    @test lbs[ind] ≈ -4
-    @test ubs[ind] ≈ 4
-
-    # Parallel test
-    wsparallel, lbsparallel, ubsparallel = warmup(
-        model,
-        Tulip.Optimizer;
-        modifications = [change_constraint("EX_glc__D_e", -4, 4)],
-        warmup_points = collect(1:n_reactions(model)),
-        workerids = W,
-    )
-
-    ind = first(indexin(["EX_glc__D_e"], reactions(model)))
-    @test size(wsparallel) == (95, 2)
-    @test size(wsparallel[1, 1]) == (95,)
-    @test lbsparallel[ind] ≈ -4
-    @test ubsparallel[ind] ≈ 4
+    idx = first(indexin(["EX_glc__D_e"], reactions(model)))
+    @test size(pts) == (95, 100)
+    @test all(pts[idx, :] .>= -2)
+    @test all(pts[idx, :] .<= 2)
 end